First passage probability assessment of stationary non-Gaussian process using the third-order polynomial transformation

Abstract

In this article, an analytical moment-based procedure is developed for estimating the first passage probability of stationary non-Gaussian structural responses for practical applications. In the procedure, an improved explicit third-order polynomial transformation (fourth-moment Gaussian transformation) is proposed, and the coefficients of the third-order polynomial transformation are first determined by the first four moments (i.e. mean, standard deviation, skewness, and kurtosis) of the structural response. The inverse transformation (the equivalent Gaussian fractile) of the third-order polynomial transformation is then used to map the marginal distributions of a non-Gaussian response into the standard Gaussian distributions. Finally, the first passage probabilities can be calculated with the consideration of the effects of clumping crossings and initial conditions. The accuracy and efficiency of the proposed transformation are demonstrated through several numerical examples for both the “softening” responses (with wider tails than Gaussian distribution; for example, kurtosis > 3) and “hardening” responses (with narrower tails; for example, kurtosis < 3). It is found that the proposed method has better accuracy for estimating the first passage probabilities than the existing methods, which provides an efficient and rational tool for the first passage probability assessment of stationary non-Gaussian process.

Keywords

first passage probability non-Gaussian structural responses the equivalent Gaussian fractile the first four moments third-order polynomial transformation

Introduction

The failure of a random dynamic system is usually defined based on the exceedance of a predefined critical threshold level. In many engineering applications, it is important to calculate the probabilities that a structural response exceeds the critical threshold levels in the given time interval [0, T]. The critical threshold level could be the yield stress, strain, or displacement. Failure occurs when the critical threshold value is exceeded at its first time, and this kind of analysis falls into the category of the first passage problem (Coleman, 1959; Lin, 1970; Rice, 1944, 1945; Vanmarcke, 1975).

So far, there are no simple analytical solutions existing for the first passage problem even for the simplest case, that is, a single-degree-of-freedom (SDOF) linear oscillator subject to Gaussian white noise excitation with a deterministic failure threshold (Crandall, 1970; Naess, 1990). Many numerical studies have been devoted to the computation of the first passage probability (Au and Beck, 2001a, 2001b; Naess et al. 2010; Naess and Gaidai, 2009), and the Monte Carlo simulation (MCS) is the most general and robust, although it is also computationally expensive for low-probability events. For analytical studies, it has also been devoted to estimate the first passage probability (Barbato and Conte, 2011; Coleman, 1959; Ghazizadeh et al., 2012; Lin, 1970; Rice, 1944, 1945; Vanmarcke, 1975), and the Poisson approximation is the most widely used. In most cases, the computation of the first passage probability is often based on the up-crossing approach which assumes that the crossings are independent and the responses are Gaussian.

The well-known Gaussian model is generally inappropriate for calculating the first passage probability of a structure when structural behavior is non-linear, or the excitation is non-Gaussian (e.g. wind and wave loads). Structural responses generally cannot be modeled as Gaussian stochastic processes for these cases. This is important particularly for non-Gaussian response whose tail behavior is different from Gaussian responses. Series distribution methods (Grigoriu, 1984; Ochi, 1986) can be used to transform known statistical results (such as mean up-crossing rates and extreme values of a Gaussian process) into those of a non-Gaussian process by finding a simple functional transformation of the equivalent Gaussian statistics. For these cases, various non-linear models have been formulated through series approximation to map a non-Gaussian response process into a standard Gaussian one. For example, full probability distributions have been estimated from response moments (e.g. mean, standard deviation, skewness, and kurtosis) with Gram–Charlier and Edgeworth series (Johnsom and Kotz, 1970). These series, however, can yield multi-modal and even negative probability densities and crossing rates. Winterstein (1988) proposed Hermite moment models known as Winterstein’s polynomials (see Appendix 1) using response moments in terms of Hermite polynomials to form non-Gaussian response contributions. Winterstein’s polynomials offers equivalent Gaussian fractile mapping between Gaussian and non-Gaussian processes and has been widely used (Jensen, 1994; Puig and Akian, 2004; Winterstein and Kashef, 2000; Winterstein and Mackenzie, 2011). Winterstein’s polynomials, however, cannot be directly applied to some processes (Choi and Sweetman, 2010) since the transformation is non-monotonic for certain combinations of skewness and kurtosis (Low, 2013) and is known to have poor accuracy for strongly non-Gaussian processes (Ding and Chen, 2014). Based on Winterstein’s polynomials, the first passage probability for stationary non-Gaussian response with the consideration of the effects of clumping and initial conditions on the up-crossing rate has been proposed by He and Zhao (2007). However, significant difference was observed between the results based on Winterstein’s polynomials and those from MCS method, especially for the case of “hardening” responses (with narrower tails; for example, kurtosis < 3) due to transformation errors. Moreover, unlike the “softening” responses (with wider than Gaussian distribution tails; for example, kurtosis > 3) of Winterstein’s polynomials, the relationship between statistical moments and Hermite model coefficients cannot be formulated in closed form (Ding and Chen, 2016). Further research is, therefore, necessary to improve the accuracy of transforming the stationary non-Gaussian processes to Gaussian ones. Although the reliability problem where the discrete load process (such as a Poisson process) is also studied in some relevant references (Wang et al., 2016; Wang and Zhang, 2018), it is should be noted that the focus of this article is on the reliability problem subject to a continuous load process.

The objective of this article is to propose a new transformation approach properly for both the “softening” and “hardening” responses to calculate the first passage probabilities of stationary non-Gaussian structural responses. This article is organized as follows. In the next section, the formulation of the first passage probabilities of stationary non-Gaussian structural responses is described; in section “The third-order polynomial transformation for the first passage probability assessment of stationary non-Gaussian process”, the third-order polynomial transformation is developed to approximate stationary non-Gaussian structural responses and for practical applications, an explicit third-order polynomial transformation (fourth-moment Gaussian transformation) is proposed based on comparison studies of transformation methods, and its inverse transformation (the equivalent Gaussian fractile) is then derived to map a non-Gaussian structural response into a standard Gaussian process; in section “Numerical examples”, the accuracy and efficiency of this transformation for both “softening” responses and “hardening” responses are demonstrated through several numerical examples; finally, the findings of this article are summarized in the conclusions. It is believed that the methodology presented in this article provides an efficient and rational tool for the first passage probability assessment of stationary non-Gaussian process.

Formulation of the first passage probability of stationary non-Gaussian structural responses

In the structural safety analysis, one often assumes that the crossing of a random response follow the Poisson distribution and the first passage probability, P_f(T) is given by (Coleman, 1959)

P_{f} (T) = 1 - \exp {- E [N^{+} (T)]}

(1)

where $E [N^{+} (T)]$ is the mean number of crossings of X(t) exceeding the critical level x during the interval [0, T] for regular streams of up-crossings. In the stationary case, one has $E [N^{+} (T)] = v^{+} T$ , where $v^{+}$ denotes the mean up-crossing rate and can be estimated by Rice (1944) formula as

v^{+} (t) = \int_{0}^{\infty} \overset{\cdot}{x} f_{X \overset{\cdot}{X}} (x, \overset{\cdot}{x}, t) d \overset{\cdot}{X}

(2)

where $f_{X \overset{\cdot}{X}} (x, \overset{\cdot}{x}, t)$ is the joint probability density function (PDF) of X(t) and its derivative $\overset{\cdot}{X} (t)$ at time t. For a stationary standard Gaussian process U(t) at critical level u, equation (2) can be written as

v^{+} = v_{0} \exp (- \frac{u^{2}}{2}), v_{0} = \frac{σ_{\overset{\cdot}{U}}}{2 π}

(3)

where $v_{0}$ is the mean rate of response cycles, $σ_{\overset{\cdot}{U}}$ is the standard deviation of $\overset{\cdot}{U} (t)$ , and $\overset{\cdot}{U} (t)$ is the derivative of U(t) at time t. For a linear oscillator, $v_{0}$ can be written as

v_{0} = \frac{ω_{0}}{2 π}

(4)

where $ω_{0}$ is the undamped natural frequency of a dynamical system.

This approach has been proved to be exact for very high critical levels (Cramér, 1966). If X(t) is narrow-banded or the critical threshold level is low, the clumping effect of the up-crossing events should be considered to modify the mean level up-crossing rates in equation (2) (He, 2009, 2015). Since the exceedances of the critical threshold level tend to occur in clumps, the clump size $〈 cs 〉$ is proposed to improve equation (1) (Langley, 1988). When the stationary start condition and the clump size are considered, equation (1) can be further improved as (Ditlevsen, 1986)

P_{f} (T) = 1 - [1 - P_{f} (0)] \exp {- \frac{v^{+} T}{〈 cs 〉 \cdot [1 - P_{f} (0)]}}

(5)

where

\begin{matrix} 〈 cs 〉 = 1 + \sum_{n = 2}^{\infty} p (cs ⩾ n), p (cs ⩾ n) \\ = 2 {1 - Φ [\frac{u (x) \cdot (1 - R_{U} (τ_{n}))}{\sqrt{1 - R_{U}^{2} (τ_{n})}}]} \end{matrix}

(6)

v^{+} = v_{0} \exp [- \frac{u^{2} (x)}{2}], P_{f} (0) = 1 - Φ [u (x)]

(7)

For high critical level, equation (1) can be improved as (Ditlevsen, 1986)

P_{f} (T) = 1 - [1 - P_{f} (0)] \exp {- \frac{v^{+} T}{[1 - P_{f} (0)]}}

(8)

where P_f(0) is the probability of instantaneous failure at T = 0; $R_{U} (τ)$ is the auto correlation function of the standard linear Gaussian process U(t), which can be determined as shown in Appendix 2, and $τ_{n} = t_{n} - t_{1}$ is a suitable choice of τ; $Φ [\cdot]$ is the cumulative distribution function (CDF) of standard normal variables; u(x) is the equivalent Gaussian fractile. Equations (5) to (8) can be considered as a general form for calculating first passage probability of a stationary non-Gaussian process.

The first passage probabilities were investigated by He and Zhao (2007) for stationary non-Gaussian structural responses using equations (5) to (8), in which the equivalent Gaussian fractile u(x) based on Winterstein’s polynomials was adopted. However, significant differences were observed between the results based on Winterstein’s polynomials and those estimated using the MCS method, especially for the case of “hardening” responses due to the transformation errors of Winterstein’s polynomials. In the next section, the third-order polynomial transformation is proposed to solve this problem.

The third-order polynomial transformation for the first passage probability assessment of stationary non-Gaussian process

The third-order polynomial transformation and its inverse function

The third-order polynomial transformation method suggested by Fleishman (1978) can be formulated as

\frac{X (t) - μ_{X}}{σ_{X}} = X_{S} (t) = a_{1} + a_{2} U (t) + a_{3} U^{2} (t) + a_{4} U^{3} (t)

(9)

where X(t) is a stochastic process; X_S(t) is standardized stochastic process; $μ_{X}$ and $σ_{X}$ are the mean value and standard deviation of X(t), respectively; and a₁, a₂, a₃, and a₄ are polynomial coefficients that can be obtained by making the first four moments of the left side of equation (9) equal to those of the right side (Fleishman, 1978; Zhao and Lu, 2008).

Once the four coefficients determined and considering the critical threshold level x, the inverse transformation (the equivalent Gaussian fractile) u(x) of the third-order polynomial transformation can be expressed as follows (Zhang et al., 2018; Zhao and Lu, 2007):

When a₄ > 0 and $p = (3 a_{2} a_{4} - a_{3}^{2}) / (3 a_{4}^{2}) ⩾ 0$

u (x) = \sqrt[3]{A} + \sqrt[3]{B} - \frac{a}{3}

(10)

where A, B, and a are the parameters, and they are expressed as follows

A = - \frac{q}{2} + \sqrt{Δ}, B = - \frac{q}{2} - \sqrt{Δ}

(11)

\begin{matrix} Δ = {(\frac{p}{3})}^{3} + {(\frac{q}{2})}^{2}, q = \frac{2}{27} a^{3} - \frac{ac}{3} - a - \frac{x_{s}}{a_{4}}, \\ a = \frac{a_{3}}{a_{4}}, c = \frac{a_{2}}{a_{4}}, x_{s} = \frac{x - μ_{X}}{σ_{X}} \end{matrix}

(12)

When a₄ > 0, $p = (3 a_{2} a_{4} - a_{3}^{2}) / (3 a_{4}^{2}) < 0$ , and $α_{3 X} ⩾ 0$

u (x) = {\begin{matrix} 2 r \cos (θ / 3) - a / 3, & J_{1}^{*} < x < J_{2}^{*} \\ \sqrt[3]{A} + \sqrt[3]{B} - a / 3, & x ⩾ J_{2}^{*} \end{matrix}

(13)

where

\begin{matrix} J_{1}^{*} = σ_{x} a_{4} (- 2 r^{3} + \frac{2}{27} a^{3} - \frac{ac}{3} - a) + μ_{x}, \\ J_{2}^{*} = σ_{x} a_{4} (2 r^{3} + \frac{2}{27} a^{3} - \frac{ac}{3} - a) + μ_{x}, r = \sqrt{- \frac{p}{3}} \end{matrix}

(14)

When a₄ > 0, $p = (3 a_{2} a_{4} - a_{3}^{2}) / (3 a_{4}^{2}) < 0$ , and $α_{3 X} < 0$

u (x) = {\begin{matrix} 2 r \cos (θ / 3) - a / 3, & J_{1}^{*} < x < J_{2}^{*} \\ \sqrt[3]{A} + \sqrt[3]{B} - a / 3, & x ⩽ J_{1}^{*} \end{matrix}

(15)

When a₄ < 0

u (x) = - 2 r \cos [(θ + π) / 3] - a / 3

(16)

When a₄ = 0 and a₃ ≠ 0

u (x) = \frac{- a_{2} + \sqrt{a_{2}^{2} + 4 a_{3} (a_{3} + x_{s})}}{2 a_{3}}

(17)

When a₄ = a₃ = 0

u (x) = x_{s}

(18)

Combined with equations (5) to (8) and the inverse transformation (the equivalent Gaussian fractile) of the third-order polynomial transformation u(x) given by equations (10), (13), and (15) to (18), an analytical moment-based procedure based on the third-order polynomial transformation is developed for calculating the first passage probability of a stationary non-Gaussian process.

An improved explicit third-order polynomial transformation and its inverse function

As the determination of the four coefficients of equation (9) is not easy, the solution of non-linear equations has to be found. Therefore, it is attractive for practical applications to have the approximate but explicit formulae for the four coefficients in equation (9).

In order to avoid the solution of non-linear equations, Zhao and Lu (2007) proposed a simple explicit third-order polynomial function, which was developed from a large amount of data of third and fourth dimensionless central moments through trial and error

X_{S} (t) = - l_{1} + k_{1} U (t) + l_{1} U^{2} (t) + k_{2} U^{3} (t)

(19)

where the coefficients $l_{1}$ , $k_{1}$ , and $k_{2}$ are given as

\begin{matrix} l_{1} = \frac{α_{3 X}}{6 (1 + 6 l_{2})}, k_{1} = \frac{1 - 3 l_{2}}{(1 + l_{1}^{2} - l_{2}^{2})}, \\ k_{2} = \frac{l_{2}}{(1 + l_{1}^{2} + 12 l_{2}^{2})} \end{matrix}

(20)

and $l_{2}$ is expressed as

l_{2} = \frac{\sqrt{6 α_{4 X} - 8 α_{3 X}^{2} - 14} - 2}{36}

(21)

It was shown that equations (19) to (21) can give fairly good approximation of equation (9) when $| {α_{3}}_{X} |$ is less or equal to 1.2 (Zhao and Lu, 2007). For most practical cases of interests (Kwon and Kareem, 2011), the ranges of skewness and kurtosis are $| α_{3 X} | < \sqrt{2 (α_{4 X} - 3) / 3}$ and 3 < α_4X < 15. When $1.2 < | α_{3 X} | < \sqrt{2 (α_{4 X} - 3) / 3}$ , it is not clear whether equations (19) to (21) can provide good approximation of equation (9) or not, therefore, requires further investigation.

The four polynomial coefficients determined by equations (19) to (21) are illustrated and compared with exact values obtained by non-linear equations (Zhao and Lu, 2007) in Figure 1. The coefficients are expressed as functions of $α_{4 X}$ for α_3X = 1.6, 2.0, and 2.4. As can be observed from Figure 1 that the coefficients a₁ and a₃ predicted by equations (20) and (21) are in close agreement with the exact ones, while the coefficients a₂ and a₄ predicted by equations (20) and (21) differ significantly from the accurate ones.

Figure 1.

Fitting and validation of polynomial coefficients for the proposed model.

In this study, an improved version of equation (20) is developed from a large amount of data of third and fourth dimensionless central moments through trial and error. The development procedure is simply described as: for a given skewness α_3X and different kurtosis α_4X of X_S(t), a set of polynomial coefficients (accurate value) can be determined by making the first four moments of the left side of equation (9) equal to those of the right side (Fleishman, 1978; Zhao and Lu, 2008); Then, using regression analysis, the polynomial coefficients as functions of skewness α_3X and kurtosis α_4X can be obtained as follows

When |α_3X| ⩽ 1.2, equations (20) and (21) are still applicable, while when $1.2 < | α_{3 X} | < \sqrt{2 (α_{4 X} - 3) / 3}$ , the coefficients l₁, k₁, and k₂ are modified as $l'_{1}$ , $k'_{1}$ , and $k'_{2}$

\begin{matrix} l'_{1} = \frac{α_{3 X}}{6 (1 + 6 l_{2} + 2 l_{2}^{2})}, k'_{1} = \frac{1 - 3 l_{2}}{(1 + 2 l_{1}'^{2} + 3 l_{2}^{2})} \end{matrix}, k'_{2} = \frac{l_{2}}{(1 - l_{1}'^{2} + 6 l_{2}^{2})}

(22)

where $l_{2}$ is given by equation (21).

The four polynomial coefficients determined by equation (22) and Winterstein’s polynomials (Winterstein, 1988) are also illustrated in Figure 1. As can be observed from Figure 1 that (1) coefficients of Winterstein’s polynomials have the largest differences from the accurate coefficients and (2) coefficients obtained using equation (22) are in close agreement with the accurate ones throughout the entire investigation range.

According to equations (10) to (18), the inverse transformation or equivalent Gaussian fractile u(x) based on the improved explicit third-order polynomial transformation (equations (19) to (22)) can be readily obtained through replacing a₁, a₂, a₃, and a₄ by −l₁, k₁, l₁, and k₂ or $- l'_{1}$ , $k'_{1}$ , $l'_{1}$ , and $k'_{2}$ , respectively.

The first passage probability assessment of stationary non-Gaussian process based on the third-order polynomial transformation

According to the description of sections “Formulation of the first passage probability of stationary non-Gaussian structural responses” and “The third-order polynomial transformation for the first passage probability assessment of stationary non-Gaussian process”, the procedure for calculating the first passage probability of a stationary non-Gaussian process based on the third-order polynomial transformation is described as follows, with the flowchart illustrated in Figure 2:

Estimate the first four moments of a stationary non-Gaussian process.

Determine the critical threshold level x of a stationary non-Gaussian process.

After the first four moments of a stationary non-Gaussian process and the critical threshold levels x have been given, the equivalent Gaussian fractile based on the explicit third-order polynomial transformation can be estimated by equations (10) to (18), where the coefficients a₁, a₂, a₃, and a₄ are replaced by −l₁, k₁, l₁, and k₂ or $- l'_{1}$ , $k'_{1}$ , $l'_{1}$ and $k'_{2}$ , which are determined by equations (20) to (22).

Compute the first passage probability of a stationary non-Gaussian process using equations (5) to (8).

Figure 2.

Flowchart of procedure for calculating the first passage probability of a stationary non-Gaussian process.

Numerical examples

To investigate the accuracy and efficiency of the proposed explicit third-order polynomial transformation, the x–u transformations of three different distributions including “softening” responses ( $| α_{3 X} | ⩽ 1.2$ and $| α_{3 X} | > 1.2$ ) and “hardening” responses (α_4X < 3) are discussed in detail in Example 1. Then, the application of the third-order polynomial transformation for the first passage probability assessment of stationary non-Gaussian process is demonstrated by a “softening” responses of a linear SDOF system, a “hardening” responses of a Duffing oscillator and a six-story moment-resisting steel frame (multi-degree-of-freedom linear system) subjected to stochastic wind load, respectively, in Examples 2, 3, and 4, in which MCS is employed to examine the precision of the proposed method.

x–u transformations using the first four moments of non-normal random variables

Example 1: the investigation of the proposed x–u transformations

In evaluating a normal transformation technique, the first concern could be how the relation between non-normal and normal variables is described by the technique. Suppose a random variable x is known to have a PDF, f(x), the x–u transformation can be obtained using the proposed fourth-moment standardization function or the other normal transformation techniques. Since the Rosenblatt transformation (Hohenbichler and Rackwitz, 1981) completely preserves the known marginal distribution, that is, F(x) = Φ(u), it is used herein as the benchmark in performance evaluation for the other normal transformation techniques.

Three random variables with different distributions, that is, Gumbel distribution, log-normal distribution, and truncated normal distribution, are considered in the first example. The parameters of the three distributions and their first four moments are listed in Table 1. From Table 1, Case 1 is a “softening” response with kurtosis = 5.4 (>3) and $| α_{3 X} | = 1.14 ⩽ 1.2$ ; Case 2 is a “softening” response with kurtosis = 12.922 (>3) and $1.2 < | α_{3 X} | = 2.225 < \sqrt{2 (α_{4 X} - 3) / 3} = 2.572$ ; and Case 3 is a “hardening” response with kurtosis = 2.794 (<3). The variations of the x–u transformation function with respect to x_s are shown in Figures 3 to 5 for Case 1, Case 2, and Case 3, respectively, for the results obtained using the Rosenblatt transformation, the proposed transformation, and Winterstein’s polynomial. Figures 3 to 5 reveal the following:

In Figure 3, the transformation function obtained using Winterstein’s polynomials provides good results when the absolute value of x_s is small, while when the absolute value of x_s is large, the results obtained from Winterstein’s polynomials differ greatly from those obtained using the Rosenblatt transformation.

In Figure 4, the transformation function based on Winterstein’s polynomials yields significant errors in almost the whole investigation range, compared with those obtained using the Rosenblatt transformation. The errors come from the transformation error of Winterstein’s polynomials especially when the absolute value of $α_{3 X}$ is large $(| α_{3 X} | = 2.225)$ as shown in Figure 1.

Again, it can be clearly observed from Figure 5 that the transformation function obtained using Winterstein’s polynomial has significant difference from those obtained using the Rosenblatt transformation when the absolute value of x_s is large.

The proposed method shows better accuracy than Winterstein’s polynomials. The results obtained from the proposed transformation coincide with those of the Rosenblatt transformation throughout the entire investigation range, which further demonstrate the accuracy and efficiency of the proposed polynomial coefficients in the previous section.

Table 1.

The probability distributions and their first four moments for Example 1.

Cases	Distributions	Parameters	µ_X	σ_X	α _3X	α _4X
Case 1	Gumbel	α = 0.550, β = 0.780	1	1	1.14	5.4
Case 2	log-normal	λ = 2.126, η = 0.594	10	6.5	2.225	12.922
Case 3	Truncated normal	µ = 2, σ = 1, a = −1, b = 3	1.717	0.785	−0.539	2.784

α and β are the shape parameter and scale parameter of Gumbel distribution, respectively; λ and η are shape parameter and scale parameter of log-normal distribution, respectively; a and b are the minimum and maximum of truncated normal distribution, respectively.

Figure 3.

x–u transformations of Gumbel distribution in Example 1.

Figure 4.

x–u transformations of log-normal distribution in Example 1.

Figure 5.

x–u transformations of truncated normal distribution in Example 1.

Application in first passage probability assessment of stationary non-Gaussian process

Example 2: first passage probabilities of softening responses

The second example considers a linear SDOF structure excited by a quadratic forcing function F(t) as shown in Figure 6. Specifically, it is assumed that the forcing function is F(t) = α₁U(t) + α₂U²(t), where U(t) is a stationary, zero-mean Gaussian process, and α₁ and α₂ are constants. This example is chosen primarily for its simplicity, but a similar forcing function could, for instance, serve as a model for dynamic wind pressure, although with parameters different from those chosen here.

Figure 6.

A single-degree-of-freedom (SDOF) structure excited by a quadratic forcing function F(t) for Example 2.

According to Naess (1987), the non-Gaussian structural response, X(t), to the excitation force F(t) is

\begin{matrix} X (t) = \int_{0}^{\infty} α_{1} h (τ) U (t - τ) d τ + \int_{0}^{\infty} α_{2} h (τ) U^{2} (t - τ) d τ \\ = \int_{0}^{\infty} α_{1} h (τ) U (t - τ) d τ + \int_{0}^{\infty} \int_{0}^{\infty} α_{2} h (τ_{1}) δ (τ_{1} - τ_{2}) \\ U (t - τ_{1}) U (t - τ_{2}) d τ_{1} d τ_{2} \\ = \int_{0}^{\infty} α_{1} h (τ) U (t - τ) d τ + \int_{0}^{\infty} \int_{0}^{\infty} h_{2} (τ_{1}, τ_{2}) \\ U (t - τ_{1}) U (t - τ_{2}) d τ_{1} d τ_{2} \end{matrix}

(23)

where $h (τ)$ is an ordinary impulse response function of linear equation of motion of the system, $δ (\cdot)$ is Dirac’s delta function, and $h_{2} (τ_{1}, τ_{2})$ is referred to as the quadratic impulse response function which characterizes the second-order properties of the physical system.

Equation (23) gives the known second-order stochastic Volterra series solution of the considered non-linear dynamic system. How to use numerical calculation to obtain statistical moments of the stationary responses of the non-linear dynamical system can be found in Naess (1987).

In this example, it is assumed that $α_{1} = α_{2} = ω_{0}^{2}$ , ω₀ = 1.0 rad/s, and the variance spectrum S_Q(ω) of the generalized load is 0.5 m²s between ω_L = 0.5 rad/s and ω_U = 1.5 rad/s, zero outside this interval for ω > 0 and S_Q(−ω) = S_Q(ω). The first four moments of quadratic part of stationary structural responses X(t), that is, the responses to α₂ < U²(t) determined under the damping ratio (ς) equals to 0.1, are µ_X = 1.0, σ_X = 1.62, α_3X = 1.7, and α_4X = 9.26, respectively (Naess, 1987) (“softening” responses, that is, α_4X > 3).

Using the first four moments of the quadratic part of the responses, the first passage probabilities of stationary non-Gaussian structural responses for critical threshold level x = 4.5σ_X and x = 5σ_X obtained using the proposed method are shown in Figure 7(a) and (b), respectively, compared with the results obtained by assuming Gaussian model, the method proposed by He and Zhao (2007), and the MCS method (Au and Beck, 2001b; Bayer and Bucher, 1999). From Figure 7(a) and (b):

The results of Gaussian model have the greatest differences from those of MCS method (with 100,000 samples) and are too conservative.

Although the method by He and Zhao (2007) improves the Gaussian model, its results are still somewhat quite different from those by MCS method. The difference comes from transformation errors of Winterstein’s polynomials, to which the first passage probability has large sensitivity.

The proposed method performs better than the method by He and Zhao (2007), and the obtained results are in close agreement with those by MCS method during the whole interval T. Therefore, for “softening” responses, the proposed procedure is more suitable for engineering applications.

Figure 7.

The first passage probabilities of structural responses in Example 2: (a) x = 4.5σ_X and (b) x = 5σ_X.

Example 3: first passage probabilities of hardening responses

The third example considers a Duffing oscillator subjected to a stationary Gaussian white noise. The structural response X(t) due to the excitation f(t) is determined by the following equation

\overset{\cdot\cdot}{X} (t) + 2 ς ω_{0} \overset{\cdot}{X} (t) + ω_{0}^{2} [X (t) + ε X^{3} (t)] = f (t)

(24)

where ε is a parameter for denoting the non-linearity of the restoring force; and f(t) is a zero-mean stationary Gaussian white noise. Using Fokker-Planck-Kolmogorov (FPK) equation method (Crandall, 1958), the PDF of non-linear displacement responses, p(X), can be expressed as

p (X) = c \exp [- \frac{1}{σ_{0}^{2}} (\frac{X^{2}}{2} + \frac{ε X^{4}}{4})]

(25)

in which

c = {\int_{- \infty}^{\infty} \exp [- \frac{1}{σ_{0}^{2}} (\frac{X^{2}}{2} + \frac{ε X^{4}}{4})] dx}^{- 1}, σ_{0}^{2} = \frac{π S_{0}}{2 ς ω_{0}^{3}}

(26)

The functional form of p(X) demonstrates that X(t) is non-Gaussian. Assuming that S₀ = 1 m²s, ω₀ = 2 rad/s, ς = 0.05, and ε = 0.5, the structural parameters and the first four moments of displacement response can be readily determined by equations (25) and (26) and are µ_X = 0, σ_X = 1.203, α_3X = 0, and α_4X = 2.369, respectively. The structural responses for this example is a “hardening” responses since α_4X = 2.369 < 3. Considering the critical threshold level of x = 2σ_X, the variations of the first passage probability with regard to excitation time T obtained by using the proposed method are illustrated in Figure 8, together with the results obtained by assuming Gaussian model, the method proposed by He and Zhao (2007), and the MCS method (Au and Beck, 2001b; Bayer and Bucher 1999) with 100,000 samples. Again, it can be clearly observed from Figure 8 that the results obtained by the proposed method are the closest to those of MCS method throughout the entire investigation range, implying that the proposed method is also applicable for the engineering cases with “hardening” responses.

Figure 8.

The first passage probabilities of non-linear single degree of freedom for x = 2σ_X in Example 3.

Example 4: first passage probabilities of a multi-degree-of-freedom linear dynamic system

The fourth example considers a six-story moment-resisting steel frame (multi-degree-of-freedom linear system) subjected to stochastic wind load as shown in Figure 9, which has been investigated by Au and Beck (2001b) and He (2010). It is should be noted that only a continuous wind load process is considered in the example, with a relatively short time period. The member sections of the frame are given in Table 2. For each of floor, the same section is used for all girders. The structure is modeled as a two-dimensional linear frame with beam elements connecting the joints of the frame. Masses including the contributions from the dead load of the floors and the frame members are lumped at the nodes of the frame and are tabulated in Table 3. The linear dynamical system is governed by the differential equation

M \overset{\cdot\cdot}{X} (t) + C \overset{\cdot}{X} (t) + KX (t) = F (t)

(27)

where X (t) is the N × 1 displacement vector, M is the N × N mass matrix, C is the N × N damping matrix, K is the N × N stiffness matrix, and F (t) is the N × 1 vector of stochastic wind load.

Figure 9.

A six-story moment-resisting frame structure for Example 4.

Table 2.

Sections for frame members for Example 4.

Story	Exteriorcolumn	Interior column	Girder
1, 2	C1: W14 × 159	C1: W27 × 161	G1: W24 × 94
3, 4	C2: W14 × 132	C2: W27 × 114	G2: W24 × 76
5, 6	C3: W14 × 99	C3: W27 × 84	G3: W24 × 55

Table 3.

Point masses for Example 4.

Floor	Exteriorcolumn (×10³ kg)	Interior column(×10³ kg)
2	60.4	80.0
3	53.3	78.1
4	51.9	76.0
5	51.7	75.8
6	50.1	73.5
Roof	44.6	63.1

The undamped natural frequencies of the first six modes are computed as 0.55166, 1.56334, 2.83893, 4.48842, 6.17028, and 8.16689 Hz. Rayleigh damping is assumed for the first two modes with 5% of critical damping.

Let the wind pressure $W = 2 c \bar{v} U (t) + c U^{2} (t)$ , where $\bar{v}$ is the mean wind speed; U(t) is a stationary standard Gaussian process; and c is a constant. For this study, it is assumed that $\bar{v} = 1 m / s$ and c = 1/1600, then $W = U (t) / 800 + U^{2} (t) / 1600$ (kN/m²); shape coefficient of building is µ_s₁ = 0.8 (windward side) and µ_s₂ = 0.5 (leeward side); the front face area A_j of each floor is 1.905 m × 7.2 m, 3.81 m × 7.2 m, 3.81 m × 7.2 m, 3.81 m × 7.2 m, 3.81 m × 7.2 m, and 3.81 m × 7.2 m (from top to bottom). Then the stochastic wind load, F(t) = α₁U(t) + α₂U²(t), is a stationary non-Gaussian stochastic process and can be written as

F (t) = (μ_{S 1} + μ_{S 2}) W A_{j}

(28)

The output non-Gaussian structural responses Y_m(t), m = 1, …, 24 denotes the inter-story drift ratios of the columns at every floor. The unit impulse response function, h_m(t), corresponding to Y_m(t) includes the dynamics of the frame structure and the soil layers, and is obtained by deterministic dynamic analysis. The first four moments of the inter-story drift ratios of the exterior columns at every floor for the total responses, that is, the response to F(t) and for the quadratic responses, that is, the responses to α₂U²(t), can be readily determined by MCS method and are listed in Table 4.

Table 4.

The first four moments of the inter-story drift ratios of the exterior columns at every floor for the total responses and for the quadratic responses for Example 4.

Floor	The first four moments for the total responses				The first four moments for the quadratic responses
Floor	µ_X	σ_X	α _3X	α _4X	µ_X	σ_X	α _3X	α _4X
2	0.42309	0.32281	1.04012	4.03776	0.29068	0.21058	0.91458	3.71255
3	0.47164	0.36548	1.08280	4.09908	0.31846	0.23554	0.94161	3.74786
4	0.48094	0.37262	1.10851	4.47168	0.31765	0.23910	1.03170	4.15948
5	0.43383	0.33991	1.17082	4.61608	0.27872	0.21128	1.05805	4.17780
6	0.40643	0.32103	1.13254	4.44469	0.24699	0.19363	1.13119	4.85518
Roof	0.13693	0.11014	1.23775	4.82547	0.08182	0.06641	1.32126	5.54932

Using the first four moments of total responses as shown in Table 4 and considering the critical threshold level of x = 5σ_X, the variations of the first passage probability of total responses of the inter-story drift ratios with regard to excitation time T obtained using the proposed method are illustrated in Figure 10(a) to (f), together with the results obtained by assuming Gaussian model, the method proposed by He and Zhao (2007), and the MCS method (Au and Beck, 2001b; Bayer and Bucher 1999) with 100,000 samples. It can be observed that results obtained by the proposed method shows better agreements with those from MCS method than those obtained by He and Zhao method and Gaussian model, which indicates that the proposed method can be used in multi-degree-of-freedom dynamical systems for practical use.

Figure 10.

The first passage probabilities of inter-story drift ratios exceeding critical threshold level x = 5σ_X for the total responses in Example 4: (a) sixth story (b) fifth story, (c) fourth story, (d) third story, (e) second story, and (f) first story.

Similarly, using the first four moments of the quadratic responses in Table 4 and considering the critical threshold level of x = 5σ_X, the variations of the first passage probability of quadratic responses of the inter-story drift ratios with regard to excitation time T obtained using the proposed method are illustrated in Figure 11(a) to (f), together with the results obtained by assuming Gaussian model, the method proposed by He and Zhao (2007), and the MCS method (Au and Beck, 2001b; Bayer and Bucher 1999) with 100,000 samples. It can also be clearly observed that the results obtained by the proposed method are closer to MCS method than those obtained by He and Zhao method and Gaussian model throughout the entire investigation range.

Figure 11.

The first passage probabilities of inter-story drift ratios exceeding critical threshold level x = 5σ_X for the quadratic responses in Example 4: (a) sixth story, (b) fifth story, (c) fourth story, (d) third story, (e) second story, and (f) first story.

Conclusion

From the investigations of this article, the following conclusions can be drawn:

A third-order polynomial transformation (fourth-moment Gaussian transformation) method is proposed for estimating the first passage probability of stationary non-Gaussian structural responses. An explicit third-order polynomial transformation is first proposed, and the inverse transformation (the equivalent Gaussian fractile) of the third-order polynomial transformation is then derived to map a non-Gaussian structural response into a standard Gaussian process.

Since the proposed explicit third-order polynomial transformation can give a good approximation for the polynomial coefficients using the first four moments, the proposed method provides more appropriate x–u transformation results compared to the current formulae.

Through numerical examples presented, it is demonstrated that for both the “softening” and “hardening” responses, the proposed method is sufficiently accurate in the calculation of first passage probabilities for stationary non-Gaussian structural responses.

The numerical results of the first passage probabilities of a moment-resisting steel frame (multi-degree-of-freedom linear dynamic system) subjected to the stochastic wind load are presented and compared with the results using other existing methods. Good agreement with MCS shows again that the proposed method is effective and accurate for practical use.

Footnotes

Appendix 1

Winterstein’s Hermite model (Winterstein, 1988) can be expressed as

(29)

\frac{X (t) - μ_{X}}{σ_{X}} = κ {U (t) + {\tilde{h}}_{3} [U^{2} (t) - 1] + {\tilde{h}}_{4} [U^{3} (t) - 3 U (t)]}

in which

(30)

\begin{matrix} {\tilde{h}}_{4} = \frac{\sqrt{1 + 1.5 (α_{4 X} - 3)} - 1}{18}, \\ {\tilde{h}}_{3} = \frac{α_{3 X}}{4 + 2 \sqrt{1 + 1.5 (α_{4 X} - 3)}} \end{matrix}

(31)

κ = \frac{1}{\sqrt{(1 + 2 {\tilde{h}}_{3}^{2} + 6 {\tilde{h}}_{4}^{2})}}

which $κ$ is a scaling factor ensuring that X(t) has unit variance; $μ_{X}$ , $σ_{X}$ , $α_{3 X}$ , and $α_{4 X}$ are the mean value, standard deviation, skewness, and kurtosis of X(t), respectively. For “softening” responses, the equivalent Gaussian fractile can be taken as

(32)

u (x) = {[\sqrt{ξ {(x)}^{2} + c} + ξ (x)]}^{1 / 3} - {[\sqrt{ξ {(x)}^{2} + c} - ξ (x)]}^{1 / 3} - a

(33)

ξ (x) = 1.5 b (a + \frac{x - μ_{X}}{κ σ_{X}}) - a^{3}

in which $a = {\tilde{h}}_{3} / 3 {\tilde{h}}_{4}^{3}$ , $b = 1 / 3 {\tilde{h}}_{4}$ , and $c = (b - 1 - a^{2})^{3}$ .

For “hardening” responses, the equivalent Gaussian fractile can be taken as

(34)

u (x) = x_{s} - \frac{α_{3 X}}{6} (x_{s}^{2} - 1) - \frac{α_{4 X} - 3}{24} (x_{s}^{3} - 3 x_{s})

Appendix 2

The mean clump size (He and Zhao, 2007) can be written as

(35)

\begin{matrix} p (cs ⩾ n) = p (X (t_{2}) > x, X (t_{3}) > x, \dots, X (t_{n}) > x | X (t_{0}) < x, X (t_{1}) > x) \end{matrix}

Using the laws of conditional probability and the fact that X(t₁ = x) implies with equal probability that X(t₁) is either the first or the last peak in a clump; if X(t) is reasonably narrow-banded, then equation (35) can be simplified as (He, 2015)

(36)

p (cs ⩾ n) = 2 p (X (t_{n}) > x | X (t_{1}) = x)

For Gaussian result, it can be shown as

(37)

p (cs ⩾ n) = 2 {1 - Φ [\frac{u (x) \cdot (1 - R_{U} (τ_{n}))}{\sqrt{1 - R_{U}^{2} (τ_{n})}}]}

where $R_{U} (τ)$ is the autocorrelation function of the standard linear Gaussian process U(t), and $τ_{n} = t_{n} - t_{1}$ is a suitable choice of τ. For linear SDOF system, $R_{U} (τ)$ can be determined by equation (38), and it can be given as (He and Zhao, 2007)

(38)

R_{U} (τ) = \exp [- ς ω_{0} τ] [\cos ω_{d} τ + (ς ω_{0} / ω_{d}) \sin ω_{d} τ]

in which $ω_{d} = ω_{0} \sqrt{1 - ς^{2}}$ , ς is damp ratio, and ω₀ is the undamped natural frequency of a dynamic system.

The probability that the clump size will be exactly equal to n is $p (cs ⩾ n) - p (cs ⩾ n + 1)$ , from which the mean clump size can be written as (He and Zhao, 2007)

(39)

\begin{matrix} 〈 cs 〉 = 1 + \sum_{n = 2}^{\infty} n [p (cs ⩾ n) - p (cs ⩾ n + 1)] \\ = 1 + \sum_{n = 2}^{\infty} p (cs ⩾ n) \end{matrix}

The infinite sum which appears in equation (39) must be truncated, as for large n. Equation (37) has the non-zero limit $2 {1 - F [u (x)]}$ . Langley (1988) suggested that the cut-off point can be calculated from

(40)

p (X (t_{n}) > x) = \exp [- u^{2} (x) / 2]

in which u(x) is the equivalent Gaussian fractile. For SDOF system, $〈 cs 〉$ is suggested by He and Zhao (2007)

(41)

〈 cs 〉 = 1 + 2 \sum_{n = 2}^{N} {1 - Φ [\frac{u (x) \cdot [1 - \exp (- 2 π ς n)]}{\sqrt{1 - \exp (- 4 π ς n)}}]}

where N is cut-off point, which is suggested as equation (42) (He and Zhao, 2007)

(42)

1 - Φ {\frac{u (x) \cdot [1 - R_{U} (τ_{n})]}{\sqrt{1 - {R_{U}}^{2} (τ_{n})}}} = \exp [- u^{2} (x) / 2]

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The research reported in this article is partially supported by the National Natural Science Foundation of China (Grant Nos 51738001, U1434204, 51422814, and 51421005), the Fundamental Research Funds for the Central Universities of Central South University (2017zzts158), and the grant of Shanghai Municipal Natural Science Foundation of China (no. 16ZR1417300). The support was gratefully acknowledged.

References

Beck

(2001a) Estimation of small failure probabilities in high dimensions by subset simulation. Probabilistic Engineering Mechanics 16(4): 263–277.

Beck

(2001b) First excursion probabilities for linear systems by very efficient importance sampling. Probabilistic Engineering Mechanics 16(3): 193–207.

Barbato

Conte

(2011) Structural reliability applications of nonstationary spectral characteristics. Journal of Engineering Mechanics 137(137): 371–382.

Bayer

Bucher

(1999) Importance sampling for first passage problems of nonlinear structures. Probabilistic Engineering Mechanics 14(1–2): 27–32.

Choi

Sweetman

(2010) The Hermite moment model for highly skewed response with application to tension leg. Journal of Offshore Mechanics and Arctic Engineering 132(2): 021602.

Coleman

(1959) Reliability of aircraft structures in resisting chance failure. Operations Research 7(5): 639–945.

Cramér

(1966) On the intersections between the trajectories of a normal stationary stochastic process and a high level. Arkiv för Matematik 6(4–5): 337–349.

Crandall

(1958) Random Vibration. Cambridge, MA: The MIT Press.

Crandall

(1970) First-crossing probabilities of the linear oscillator. Journal of Sound and Vibration 12(3): 285–299.

10.

Ding

Chen

(2014) Assessment of methods for extreme value analysis of non-Gaussian wind effects with short-term time history samples. Engineering Structures 80: 75–88.

11.

Ding

Chen

(2016) Moment-based translation model for hardening non-Gaussian response processes. Journal of Engineering Mechanics 142(2): 06015006.

12.

Ditlevsen

(1986) Duration of visit to critical set by Gaussian process. Probabilistic Engineering Mechanics 1(2): 82–93.

13.

Fleishman

(1978) A method for simulating non-normal distributions. Psychometrika 4(43): 521–532.

14.

Ghazizadeh

Barbato

Tubaldi

(2012) New analytical solution of the first passage reliability problem for linear oscillator. Journal of Engineering Mechanics 138(6): 695–706.

15.

Grigoriu

(1984) Crossings of non-Gaussian translation processes. Journal of Engineering Mechanics 110(4): 610–620.

16.

(2009) Numerical calculation for first excursion probabilities of linear systems. Probabilistic Engineering Mechanics 24(3): 418–425.

17.

(2010) An efficient numerical method for estimating reliabilities of linear structures under fully nonstationary earthquake. Structural Safety 32(3): 200–208.

18.

(2015) Approximate method for estimating extreme value responses of nonlinear stochastic dynamic systems. Journal of Engineering Mechanics 141(7): 04015009.

19.

Zhao

(2007) First passage times of stationary non-Gaussian structural responses. Computers & Structures 85(7): 431–436.

20.

Hohenbichler

Rackwitz

(1981) Zon-normal dependent vectors in structural safety. Journal of the Engineering Mechanics Division 107(6): 1227–1238.

21.

Jensen

(1994) Dynamic amplification of offshore steel platform responses due to non-Gaussian wave loads. Marine Structures 7(1): 91–105.

22.

Johnsom

Kotz

(1970) Continuous Univariate Distribution-1. Boston, MA: Houghton Mifflin Company.

23.

Kwon

Kareem

(2011) Peak factors for non-Gaussian load effects revisited. Journal of Structural Engineering 137(12): 1611–1619.

24.

Langley

(1988) A first passage approximation for normal stationary random processes. Journal of Sound and Vibration 122(2): 261–275.

25.

Lin

(1970) First-excursion failure of randomly excited structures. AIAA Journal 8(4): 720–725.

26.

Low

(2013)) A new distribution for fitting four moments and its applications to reliability analysis. Structural Safety 42(3): 12–25.

27.

Naess

(1987) The response statistics of nonlinear, second order transformations to Gaussian loads. Journal of Sound and Vibration 115(1): 103–129.

28.

Naess

(1990) Approximation first passage and extremes of narrow-band Gaussian and non-Gaussian random vibrations. Journal of Sound and Vibration 138(3): 365–380.

29.

Naess

Gaidai

(2009) Estimation of extreme values from sampled time series. Structural Safety 31(4): 325–334.

30.

Naess

Gaidai

Batsevych

(2010) Prediction of extreme response statistic of narrow-band vibrations. Journal of Engineering Mechanics 136(3): 290–298.

31.

Ochi

(1986) Non-Gaussian random processes in ocean engineering. Probabilistic Engineering Mechanics 1(1): 28–39.

32.

Puig

Akian

(2004) Non-Gaussian simulation using Hermite polynomials expansion and maximum entropy principle. Probabilistic Engineering Mechanics 19(4): 293–305.

33.

Rice

(1944) Mathematical analysis of random noise. Bell System Technical Journal 23(3): 282–332.

34.

Rice

(1945) Mathematical analysis of random noise. part III: statistical properties of random noise currents. Bell System Technical Journal 24(1): 46–156.

35.

Vanmarcke

(1975) On the distribution of the first passage time for normal stationary process. Journal of Applied Mechanics 42: 215–220.

36.

Wang

Zhang

(2018) Roles of load temporal correlation and deterioration-load dependency in structural time-dependent reliability. Computers & Structures 194: 48–59.

37.

Wang

Ellingwood

(2016) Time-dependent reliability of aging structures: an approximate approach. Structure and Infrastructure Engineering 12(12): 1566–1572.

38.

Winterstein

(1988) Nonlinear vibration models for extremes and fatigue. Journal of Engineering Mechanics 114(10): 1772–1790.

39.

Winterstein

Kashef

(2000) Moment-based load and response models with wind engineering applications. Journal of Solar Energy Engineering 122(3): 122–128.

40.

Winterstein

Mackenzie

(2011) Extremes of nonlinear vibration: comparing models based on moments, L-moments, and maximum entropy. Journal of Offshore Mechanics and Arctic Engineering 135(2): 185–195.

41.

Zhang

Zhao

(2018) The inverse transformation of the explicit fourth-moment standardization for structural reliability. Advances in Structural Engineering 21(5): 769–782.

42.

Zhao

(2007) Fourth-moment standardization for structural reliability assessment. Journal of Structural Engineering 133(7): 916–924.

43.

Zhao

(2008) Cubic normal distribution and its significance in structural reliability. Structural Engineering and Mechanics 28(3): 263–280.