L-Moments-Based FORM Method for Structural Reliability Analysis Considering Correlated Input Random Variables

Abstract

Leveraging the properties of the first three linear moments (L-moments), this study proposes an effective normal transformation for structural reliability analysis considering correlated input random variables, in which the admissible range of the initial correlation matrix when employing this transformation is also presented. Subsequently, a practical procedure for structural reliability analysis, grounded in the proposed transformation and first-order reliability method (FROM), is proposed, accommodating instances wherein the joint probability density functions (PDFs) or marginal PDFs of the relevant random variables remain unknown. In comparison to the technique premised on the first three central moments (C-moments), the proposed method, based on the first three L-moments, exhibits a more extensive applicability. Various practical scenarios showcase the method’s effectiveness and precision in calculating the structural reliability index, considering diverse distributions, numerous variables, and complex structures.

1. Introduction

As a technique for quantitatively measuring the risks corresponding to the failure events, reliability analysis plays an important role in probabilistic engineering mechanics. The quantification of structural failure probability presents a formidable challenge, leading to the development of various methods [1,2,3,4,5,6], such as first-order reliability method (FORM) [7], second-order reliability method (SORM) [8], moment-based approaches [9], and simulation methods [3].

Among these, the first-order reliability method (FORM) [7,8,9,10] has emerged as the most prevalent technique. Nevertheless, engineering practices frequently entail non-normally distributed and correlated random variables [11], thus demanding a manageable method for converting these into independent standard normal variables. In FORM, the Rosenblatt transformation [6,11] or Nataf transformation [12,13] are commonly employed for the normal transformation, provided that the random variables exhibit known CDFs or PDFs. However, in real-world applications, joint CDFs/PDFs and marginal PDFs of the probability distributions often remain unknown due to data scarcity. Moreover, the available information is always expressed as statistical data or statistical moments or as a correlation matrix.

Previous research on structural reliability with unknown probability distribution primarily concentrates on the second-order method [1,7,14,15]. Subsequently, Zhao [16] introduced a third-moment standardization function, utilizing the first three central moments (C-moments) to address independent random variables. Lu [17] subsequently enhanced the third-moment transformation to accommodate correlated random variables. However, this approach heavily depends on the first three C-moments and may exhibit inaccuracies with small sample sizes [18]. Additionally, the skewness of each original random variable is confined to a specific range [16], potentially neglecting other probability distributions.

Linear moments (L-moments), expressed as linear combinations of order statistics [18], exhibit greater robustness than C-moments, as their estimates are linear combinations of order statistics for all orders, whereas C-moments estimate amplify the powers of deviations from the mean as the order increases [19]. L-moments have demonstrated considerable improvement over C-moments, particularly for higher skewness [20,21], and have been applied across various research domains, including the polynomial normal transform of L-moments [22,23]. However, the normal transformation based on the first three L-moments has not been thoroughly examined as an efficacious instrument in structural reliability analysis. In addition, the explicit expression of equivalent correlation coefficients based on a normal transformation using the first L-moments has not been proposed.

Therefore, this study aims to develop a normal transformation method based on the first three L-moments and proposes a novel L-moments-based FORM method for structural reliability analysis considering correlated input random variables. Additionally, the application scope of the proposed normal transformation based on the first three L-moments has been investigated and discussed. The remainder of this paper is organized as follows. The first three L-moments of some common probability distributions are reviewed in Section 2. L-moments-based Normal Transformation involving the correlated variables is derived in Section 3. In Section 4, the proposed method for reliability analysis combined L-moments with FORM is proposed. Then, the admissible range of the correlation coefficients obtained by the normal transformation based on third L-moments is discussed in Section 5. Section 6, several numerical examples involving correlated variables are presented to demonstrate the accuracy and efficiency of the reliability index calculated by the proposed method, in which the results have been compared with those from MCS [24,25]. Finally, Section 7 summarizes the conclusions of the paper.

2. Review of First Three L-Moments of Some Common Probability Distributions

Let X represent a random variable characterized by a cumulative distribution function (CDF) F(x) and a quantile function x(F). Assume a random sample of size n is extracted from the distribution of X, with X_l:n < X_2:n < … < X_n:n denoting the order statistics of the sample. The quantities delineated in [18] can be employed to define the L-moments of X:

$${\lambda }_s=\frac{1}{s}{\displaystyle \sum _{k=0}^{s-1}\left[\left(-1\right)E\left(X_{s-k:s}\right)\left(\begin{array}{c}s-1\\ k\end{array}\right)\right]}\begin{array}{cc}{}^k,& s\ge 1\end{array}$$

(1)

Here, E(X_s−k:s) denotes the expectation of the (s−k)th order statistics from a set of s observations, and it can be expressed as:

$$E\left(X_{s-k}{{}_{:}}_s\right)=\frac{s!}{k!\left(s-k\right)!}{\displaystyle \int xf(x){\left[F\left(x\right)\right]}^k{\left[1-F\left(x\right)\right]}^{s-k-1}}dx$$

(2)

Substituting Equation (2) into Equation (1), an expression for the first three L-moments can be obtained as follows:

$${\lambda }_1=E\left(X_{1:1}\right)={\displaystyle {\int }_0^{1}x}dF\left(x\right)$$

(3)

$${\lambda }_2=\frac{1}{2}E\left(X_{2:2}-X_{1:2}\right)={\displaystyle {\int }_0^1\left[2F\left(x\right)-1\right]x}dF\left(x\right)$$

(4)

$${\lambda }_3=\frac{1}{3}E\left(X_{3:3}+X_{1:3}-2X_{2:3}\right)={\displaystyle {\int }_0^1\left[6F^2\left(x\right)-6F\left(x\right)+1\right]x}dF\left(x\right)$$

(5)

The concept of probability-weighted moments (PWM) of a random variable X firstly introduced by [1] can be defined as:

$$M_{p,m,r}=E\left\{X^p{\left[1-F\left(X\right)\right]}^r{\left[F\left(X\right)\right]}^m\right\}={\displaystyle {\int }_0^1{x}^p}{\left[1-F\left(x\right)\right]}^r{F}^m\left(x\right)dF\left(x\right)$$

(6)

Given a particular type of the PWM, assuming β = M_1,m,0, the following expression can be derived:

$${\beta }_m={\displaystyle {\int }_0^{1}x{F}^m\left(x\right)dF\left(x\right)}$$

(7)

It can be observed that the L-moments are essentially linear combinations of PWMs. From Equation (7) and Equations (3)–(5), the first three L-moments can be obtained as:

$${\lambda }_1={\beta }_0$$

(8)

$${\lambda }_2=2{\beta }_1-{\beta }_0$$

(9)

$${\lambda }_3=6{\beta }_2-6{\beta }_1+{\beta }_0$$

(10)

Similar to the C-moments, the normalized λ₃ can be expressed as:

$${\tau }_3={\lambda }_3/{\lambda }_2$$

(11)

where λ₁, λ₂, and τ₃ are considered to be measures of location, scale, and skewness of a probability distribution [18]. These parameters are also commonly referred to as the mean or L-location, L-scale, and L-skewness, respectively.

Table 1. The first three L-moments of some common distributions in structural reliability.

Distribution	CDFs (F(x))	L-Location (λ1)	L-Scale (λ2)	L-Skewness (τ3)
Uniform	$F\left(x\right)=\frac{x-\alpha }{\beta -\alpha }$	$\frac{\alpha +\beta }{2}$	$\frac{\beta -\alpha }{6}$	$0$
Exponential	$F\left(x\right)=1-\exp \left(-x\lambda \right)$	$\frac{1}{\lambda }$	$\frac{1}{2\lambda }$	$\frac{1}{3}$
Normal	$F\left(x\right)=\mathsf{\Phi }\left(\frac{x-\mu }{\sigma }\right)$	$\mu$	$\frac{\sigma }{\sqrt{\pi }}$	$0$
Lognormal	$F\left(x\right)=\mathsf{\Phi }\left[\frac{-\lambda +\ln \left(x-\xi \right)}{\zeta }\right]$	$\xi +\exp \left(\lambda +\frac{{\zeta }^2}{2}\right)$	$erf\left(\frac{\zeta }{2}\right)\exp \left(\lambda +\frac{{\zeta }^2}{2}\right)$	$\frac{6}{\sqrt{\pi }erf\left(\frac{\zeta }{2}\right)}{\displaystyle {\int }_0^{\frac{\zeta }{2}}erf\left(\frac{x}{\sqrt{3}}\right)}\exp \left(-x^2\right)dx$
Gamma	$F\left(x\right)=\gamma \left(\alpha ,\frac{x}{\beta }\right)\frac{1}{\Gamma \left(\alpha \right)}$	$\alpha \beta$	$\beta \frac{\Gamma \left(\alpha +\frac{1}{2}\right)}{\sqrt{\pi }\Gamma (\alpha )}$	$6I_{1/3}\left(\alpha ,2\alpha \right)-3$
Gumbel	$F\left(x\right)=\exp \left[-\exp \left(\frac{-x+\xi }{\theta }\right)\right]$	$\xi +0.5772\theta$	$0.69315\theta$	$0.1699$

provides a list of the first three L-moments for the commonly used distributions in structural reliability analysis, for ease of reference.

When estimating L-moments in practical applications, a random sample is often drawn from an unknown distribution. Let x_l, x₂, …, x_n denote the sample, and x_l:n ≤ x_2:n ≤ … ≤ x_n:n denote the ordered sample. The first three sample L-moments, i.e., l₁, l₂, and l₃, can be expressed as follows:

$$l_1=b_0$$

(12)

$$l_2=2b_1-b_0$$

(13)

$$l_3=6b_2-6b_1+b_0$$

(14)

where b₀, b₁, and b₂ are the sample PWMs corresponding to β₀, β₁, and β₂, respectively, and they can be estimated using the following expressions [19]

$$b_0=\frac{1}{n}{\displaystyle \sum _{i=1}^n{x}_{i:n}}, b_1=\frac{1}{n}{\displaystyle \sum _{i=2}^n\frac{\left(i-1\right)}{\left(n-1\right)}{x}_{i:n}}, b_2=\frac{1}{n}{\displaystyle \sum _{i=3}^n\frac{\left(i-1\right)\left(i-2\right)}{\left(n-1\right)\left(n-2\right)}{x}_{i:n}}$$

(15)

3. L-Moments-Based Normal Transformation

Owing to the incomplete probability information of random variables in practical engineering, access to limited statistical data is common. Therefore, in this section, the paper pioneeringly introduces comprehensive normal transformation and inverse transformation formulas based on the first three L-moments, along with the derivation of formulas accounting for correlated variables. These transformation formulas enable the calculation of accurate polynomial coefficients, allowing random variables to be approximated based on polynomials containing standard normal random variables. Consequently, the simulation of their actual probability distributions is made possible, which plays a significant role in the subsequent computation of structural reliability.

3.1. U-X Transformation Technique

According to the transformation technique proposed by [16], it is possible to approximate a univariate random variable X through a second-order polynomial of standard normal random variables. The approximation can be expressed as:

$$X=S_U\left(U\right)=a+bU+c{U}^2$$

(16)

The expression relates to a univariate random variable X and a standard normal random variable U. The coefficients a, b, and c are determined such that the first three L-moments of a polynomial function $S_U\left(U\right)$ match those of X.

Suppose $F(x)=\mathsf{\Phi }(u)$, and substituting Equation (16) in Equation (7), one obtains

$$\begin{array}{ll}{\beta }_m\left[S_U\left(U\right)\right]& ={\displaystyle {\int }_0^{1}x{F}^m\left(x\right)dF\left(x\right)}\\ & ={\displaystyle {\int }_{-\infty }^{\infty }{\mathsf{\Phi }}^m\left(u\right)\varphi \left(u\right)\left(a+bu+c{u}^2\right)du}, m=0, 1, 2\\ & =a{C}_{m,0}+b{C}_{m,1}+c{C}_{m,2}\end{array}$$

(17)

where $\mathsf{\Phi }(\cdot )$ and $\varphi (\cdot )$ respectively represent the CDFs and PDFs of the standard normal variable; m takes values 0, 1, or 2; and C_m,n is given by follows [26]:

$$C_{m,n}={\displaystyle {\int }_{-\infty }^{\infty }{u}^n\varphi \left(u\right){\mathsf{\Phi }}^m\left(u\right)du}, m, n=0, 1, 2$$

(18)

Thus

$$C_{0,0}=1, C_{0,1}=0, C_{0,2}=1$$

(19)

$$C_{1,0}=\frac{1}{2}, C_{1,1}=\frac{1}{2\sqrt{\pi }}, C_{1,2}=\frac{1}{2}$$

(20)

$$C_{2,0}=\frac{1}{3}, C_{2,1}=\frac{1}{2\sqrt{\pi }}, C_{2,2}=\frac{\sqrt{3}+2\pi }{6\pi }$$

(21)

Based on Equation (17) and Equations (8)–(10), the initial three L-moments of $S_U\left(U\right)$ can be computed as follows:

$${\lambda }_{1S}={\beta }_0\left[S_U(U)\right]=a+c$$

(22)

$${\lambda }_{2S}=2{\beta }_1\left[S_U\left(U\right)\right]-{\beta }_0\left[S_U\left(U\right)\right]=\frac{b}{\sqrt{\pi }}$$

(23)

$${\lambda }_{3S}=6{\beta }_2\left[S_U\left(U\right)\right]-6{\beta }_1\left[S_U\left(U\right)\right]+{\beta }_0\left[S_U\left(U\right)\right]=\frac{\sqrt{3}c}{\pi }$$

(24)

Therefore, explicit and straightforward expressions for the polynomial coefficients a, b, and c can be derived as:

$$a={\lambda }_{1X}-\frac{\pi }{\sqrt{3}}{\lambda }_{3X}$$

(25)

$$b=\sqrt{\pi }{\lambda }_{2X}$$

(26)

$$c=\frac{\pi }{\sqrt{3}}{\lambda }_{3X}$$

(27)

As indicated in Equation (16) and Equations (25)–(27), the transformation technique based on the first three L-moments does not impose any restrictions on the L-skewness. Specifically, when ${\lambda }_{1X}={\lambda }_{3X}=0$ and ${\lambda }_{2X}=\frac{1}{\sqrt{\pi }}$, the polynomial coefficients a and c both equal 0, while b equals 1, thereby resulting in the degeneration of the u-x transformation function to x = u.

3.2. X-U Transformation Technique

Based on Equation (16), the x-u transformation can be expressed as follows:

(1)

${\lambda }_{3X}\ne 0$,

$$U=\frac{-b+\sqrt{b^2-4c(a-X)}}{2c}=\frac{-\sqrt{3}{\lambda }_{2X}+\sqrt{3{\lambda }_{2X}^2+4\sqrt{3}{\lambda }_{3X}(X-{\lambda }_{1X})+4\pi {\lambda }_{3X}^2}}{2{\lambda }_{3X}}$$

(28)

(2)

${\lambda }_{3X}=0, \mathrm{and} {\lambda }_{2X}\ne 0$,

$$U=\frac{X-a}{b}=\frac{X-{\lambda }_{1X}}{\sqrt{\pi }{\lambda }_{2X}}$$

(29)

3.3. The Equivalent Correlation Coefficient

X_i and X_j are two correlated non-normal random variables. Their first three L-moments and the correlation coefficient (ρ_ij) can be used to approximate them as follows:

$$X_i=a_i+b_i{Z}_i+c_i{Z}_i^2={\left(1,Z_i,Z_i^2\right)}^T\cdot \left(a_i,b_i,c_i\right)$$

(30)

$$X_j=a_j+b_j{Z}_j+c_j{Z}_j^2={\left(1,Z_j,Z_j^2\right)}^T\cdot \left(a_j,b_j,c_j\right)$$

(31)

where Z_i and Z_j represent two correlated standard normal variables with a correlation coefficient of ρ_0ij. Then, the coefficients a_i, b_i, and c_i in Equation (30), and a_j, b_j, and c_j in Equation (31) can be obtained using Equations (25)–(27) as follows:

$$a_i={\lambda }_{1X_i}-\frac{\pi }{\sqrt{3}}{\lambda }_{3X_i}, b_i=\sqrt{\pi }{\lambda }_{2X_i}, c_i=\frac{\pi }{\sqrt{3}}{\lambda }_{3X_i}$$

(32)

$$a_j={\lambda }_{1X_j}-\frac{\pi }{\sqrt{3}}{\lambda }_{3X_j}, b_j=\sqrt{\pi }{\lambda }_{2X_j}, c_j=\frac{\pi }{\sqrt{3}}{\lambda }_{3X_j}$$

(33)

The first three L-moments of X_i and X_j, can be denoted as ${\lambda }_{1X_i}$, ${\lambda }_{2X_i}$, and ${\lambda }_{3X_i}$ and ${\lambda }_{1X_j}$, ${\lambda }_{2X_j}$, and ${\lambda }_{3X_j}$, respectively.

The correlation coefficient between Z_i and Z_j, i.e., denoted as ρ_0ij, can be obtained from the definition of the correlation coefficient. This can be expressed as follows:

$${\rho }_{ij}=\frac{Cov\left(X_i{X}_j\right)}{\sqrt{D\left(X_i\right)D\left(X_j\right)}}=\frac{E\left(X_i{X}_j\right)-E\left(X_i\right)E\left(X_j\right)}{{\sigma }_{X_i}{\sigma }_{X_j}}=\frac{E\left(X_i{X}_j\right)-{\lambda }_{1X_i}{\lambda }_{2X_j}}{{\sigma }_{X_i}{\sigma }_{X_j}}$$

(34)

where $D\left(X_i\right)$, $D\left(X_j\right)$ and ${\sigma }_{X_i}$, ${\sigma }_{X_j}$ respectively represent the variance and the standard deviation of X_i and X_j; $E\left(X_i{X}_j\right)$ can be computed as follows:

$$\begin{array}{ll}E\left(X_i{X}_j\right)& =E\left[\left(a_i,b_i,c_i\right){\left(1,Z_i,Z_i^2\right)}^T\cdot {\left(a_j,b_j,c_j\right)}^T\left(1,Z_j,Z_j^2\right)\right]\\ & =\left(a_i,b_i,c_i\right)\left[\begin{array}{ccc}1& 0& 1\\ 0& {\rho }_{0ij}& 0\\ 1& 0& 2{\rho }_{0ij}^2+1\end{array}\right]{\left(a_j,b_j,c_j\right)}^T\end{array}$$

(35)

Substituting Equations (32), (33) and (35) into Equation (34), the correlation coefficient between X_i and X_j, denoted as ρ_ij, can be expressed as a function of the correlation coefficient between the two correlated standard normal variables, denoted as ρ_0ij.

$${\rho }_{ij}=\frac{2c_i{c}_j{\rho }_{0ij}^2+b_i{b}_j{\rho }_{0ij}}{{\sigma }_{X_i}{\sigma }_{X_j}}$$

(36)

that is

$${\rho }_{ij}=\frac{2{\pi }^2{\lambda }_{3X_i}{\lambda }_{3X_j}{\rho }_{0ij}^2+3\pi {\lambda }_{2X_i}{\lambda }_{2X_j}{\rho }_{0ij}}{3{\sigma }_{X_i}{\sigma }_{X_j}}$$

(37)

The correlation coefficient, ρ_0ij, can be obtained by solving Equation (37). However, it is important to note that the obtained value of ρ_0ij must satisfy the following conditions to meet the definition of the correlation coefficient:

$$-1\le {\rho }_{ij}\le 1,-1\le {\rho }_{0ij}\le 1, {\rho }_{ij}\cdot {\rho }_{0ij}\ge 0, \left|{\rho }_{0ij}\right|\ge \left|{\rho }_{ij}\right|, \mathrm{and} \frac{\partial {\rho }_{ij}}{\partial {\rho }_{0ij}}\ge 0$$

(38)

By utilizing Equations (37) and (38),

Table 2. The formulae of ρ_0ij for different cases.

Sign of A	ρ0ij	The Range of A and B	The Upper and Lower Bounds of ρij
$A>0$	${\rho }_{0ij}=\frac{\sqrt{B^2+4A{\rho }_{ij}}-B}{2A}$	$B^2-4A\ge 0$	${\rho }_{ij-\max }=\min \left\{1,A+B\right\}$,${\rho }_{ij-\min }=\max \left\{-1,A-B\right\}$
		$B^2-4A<0$	${\rho }_{ij-\max }=\min \left\{1,A+B\right\}$,${\rho }_{ij-\min }=\max \left\{\frac{B^2}{-4A},A-B\right\}$
$A<0$	${\rho }_{0ij}=\frac{\sqrt{B^2+4A{\rho }_{ij}}-B}{2A}$	$B^2+4A\ge 0$	${\rho }_{ij-\max }=\min \left\{1,A+B\right\}$,${\rho }_{ij-\min }=\max \left\{-1,A-B\right\}$
		$B^2+4A<0$	${\rho }_{ij-\max }=\min \left\{\frac{B^2}{-4A},A+B\right\}$,${\rho }_{ij-\min }=\max \left\{-1,A-B\right\}$
$A=0$	${\rho }_{0ij}=\frac{{\rho }_{ij}}{B}$	-------	${\rho }_{ij-\max }=\min \left\{1,B\right\}$,${\rho }_{ij-\min }=\max \left\{-1,-B\right\}$

In which the operator $A=2c_i{c}_j/{\sigma }_i{\sigma }_j$ and $B=b_i{b}_j/{\sigma }_i{\sigma }_j>0$.

summarizes the expressions of the correlation coefficient ρ_0ij and the upper and lower bounds of the original correlation coefficient ρ_ij that ensure the transformation is executable. Example 1 provides an illustration of the derivation of Table 2.

Each ρ_ij-max and ρ_ij-min specify the lower and upper bounds, respectively, of the initial correlation matrix, represented as ρ_max and ρ_min. The boundary Θ of the initial correlation matrix ρ, which ensures that every element of the matrix can be applied to the proposed method, can be expressed as follows:

$$\Theta =\left[\max \left({\rho }_{ij-\min }\right),\min \left({\rho }_{ij-\max }\right)\right]$$

(39)

3.4. U-X Transformation Considering the Correlated Variables

The polynomial coefficients for the n-dimensional case can be determined using Equations (25)–(27). For any two correlated non-normal random variables in the vector, the corresponding equivalent correlation coefficients of standard normal variables can be obtained from Table 2, subject to the boundary conditions of the initial correlation matrix. This approach provides a practical and efficient technique for transforming correlated non-normal random variables to standard normal variables, which can be useful in various fields of engineering and science. The correlation matrix of the corresponding standard normal variables, C_Z, can be obtained as follows:

$$C_{\mathbf{Z}}=\left[\begin{array}{ccccc}1& \cdots & {\rho }_{01j}& \cdots & {\rho }_{01n}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\rho }_{0i1}& \cdots & {\rho }_{0ij}& \cdots & {\rho }_{0in}\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ {\rho }_{0n1}& \cdots & {\rho }_{0nj}& \cdots & 1\end{array}\right]$$

(40)

Using the Cholesky decomposition [2], the correlated standard normal random vector Z can be converted into an independent standard normal vector U, i.e., which is mathematically represented as:

$$Z=L_0\mathbf{U}$$

(41)

The lower triangular matrix, L₀, obtained from the Cholesky decomposition can be described as follows:

$$L_0=\left[\begin{array}{ccccc}{l}_{11}& \cdots & 0& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ l_{i1}& \cdots & l_{ii}& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ l_{n1}& \cdots & l_{ni}& \cdots & l_{nn}\end{array}\right]$$

(42)

In theory, the equivalent correlation matrix (C_Z) should be positive semi-definite after excluding fully correlated variables. However, computational errors during the transformation from correlated non-normal random vectors to correlated normal ones, especially for highly non-normal random variables, may result in small negative eigenvalues of C_Z. To address this issue, a method proposed by [27] is adopted. In such cases, C_Z can be expressed as follows:

$$C_Z=\mathbf{V}\mathsf{\Lambda }{\mathbf{V}}^T$$

(43)

The eigenvector matrix, V, and diagonal eigenvalue matrix, Λ, of C_Z are given by the decomposition. In cases where Λ contains small negative eigenvalues, which may arise due to computational errors during the transformation from correlated non-normal random vectors to correlated normal ones, these eigenvalues are replaced with small positive values such as 0.001. This guarantees that the Cholesky decomposition is carried out correctly.

Based on Equations (41) and (42), Z_i can be represented as follows:

$$Z_i={\displaystyle \sum _{k=1}^i{U}_k{l}_{ik}} , i=1, 2, \cdots , n$$

(44)

Upon substitution of Equation (44) into Equation (30), the transformation of correlated random variables from the U to X domain can be achieved through a second-order normal polynomial transformation, and can be expressed as follows:

$$X_i=a_i+b_i{\displaystyle \sum _{k=1}^i{U}_k{l}_{ik}}+c_i{\left({\displaystyle \sum _{k=1}^i{U}_k{l}_{ik}}\right)}^2,i=1, 2, \cdots , n$$

(45)

3.5. X-U Transformation Considering the Correlated Variables

The acquisition of the normal transformation involves first transforming the correlated non-normal random vector X into a correlated standard normal random vector Z and then transforming Z into an independent standard normal random vector U.

The X-Z transformation can be expressed using Equations (30) and (31) as follows:

$$Z_i=S_Z^{-1}\left(X_i\right)=\{\begin{array}{l}\frac{-\sqrt{3}{\lambda }_{2X_i}+\sqrt{3{\lambda }_{2X_i}^2+4\sqrt{3}{\lambda }_{3X_i}\left(X_i-{\lambda }_{1X_i}\right)+4\pi {\lambda }_{3X_i}^2}}{2{\lambda }_{3X}}, c_i\ne 0\\ \frac{X_i-{\lambda }_{1X_i}}{\sqrt{\pi }{\lambda }_{2X_i}}, c_i=0\end{array}$$

(46)

After the X-Z transformation, the correlated vector Z can be further transformed into an independent standard normal vector U,

$$U=Z{L}_0^{-1}$$

(47)

The inverse matrix L₀⁻¹ of L₀, also a lower triangular matrix, can be expressed as:

$$L_0^{-1}=\left[\begin{array}{ccccc}{h}_{11}& \cdots & 0& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ h_{i1}& \cdots & h_{ii}& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ h_{n1}& \cdots & h_{ni}& \cdots & h_{nn}\end{array}\right]$$

(48)

According to Equations (47) and (48), U_i is expressed as:

$$U_i={\displaystyle \sum _{k=1}^i{h}_{ik}{Z}_k} , i=1, 2, \cdots , n$$

(49)

The X-U transformation can be expressed as Equation (49) by substituting Equation (46) into it:

$$U_i={\displaystyle \sum _{k=1}^i{h}_{ik}{S}_Z^{-1}\left(X_k\right)}, i=1, 2, \cdots , n$$

(50)

4. The Proposed Method for Reliability Analysis Combined L-Moments with FORM

Building upon the aforementioned work, the random variables in structural analysis can be approximated as a polynomial form based on standard normal variables using the formulas proposed in the third section, even when their probability distributions remain uncertain. Consequently, in the subsequent structural reliability analysis, the situation with unknown variable distributions is transformed into a common case with known distributions (i.e., standard normal space).

In this Section, the procedure for conducting reliability analysis of a performance function that involves correlated random variables can be described using the first three L-moments, standard deviations, and correlation matrix through the FORM method. The specific steps are shown below:

(1)

Obtain the initial correlation matrix $C_X$, standard deviation, and L-moments of the correlated random variables from the correlation matrix, PDFs, or statistical data.

(2)

Using Equation (30) and Table 2, determine the polynomial coefficients in Equation (32) and the correlation coefficients between Z_i and Z_j. Then, calculate the inverse matrix L₀⁻¹ and its lower triangular matrix L₀.

(3)

Select the initial checking point x₀, which is typically the mean vector of X.

(4)

Obtain the initial reliability index β₀ by obtaining the corresponding checking point in the independent standard normal space, u₀, through the use of Equation (50).

(5)

Calculate the Jacobian matrix $J=\frac{\partial X}{\partial U}$ at the point u₀ using the elements derived from Equation (46).

$$\frac{\partial x_i}{\partial u_j}=\begin{array}{cc}{l}_{ij}\left(b_i+2c_i{\displaystyle \sum _{k=1}^i{u}_k{l}_{ik}}\right),& i,j=1,2,\dots ,n\end{array}$$

(51)

(6)

Calculate the performance function and its gradient vector at the point u₀:

$$\begin{array}{cc}g\left(u_0\right)=G\left(x_0\right),& \nabla g\left(u_0\right)=\nabla G\left(x_0\right)\cdot J^T\end{array}$$

(52)

where $\nabla g\left(u_0\right)$ represents the gradient vector of $g\left(u\right)$ evaluated at u₀, and $\nabla G\left(x_0\right)$ represents the gradient vector of G(x) evaluated at x₀.

(7)

Calculate the next design point u₁ based on Equation (53) in the independent standard normal space:

$$u_1=\frac{1}{\nabla g(u_0)\cdot {\nabla }^{T}g(u_0)}\left[u_0\cdot {\nabla }^{T}g(u_0)-g(u_0)\right]\cdot \nabla g(u_0)$$

(53)

(8)

Determine the corresponding reliability index based on Equation (54):

$$\beta =\sqrt{u_1^T{u}_1}$$

(54)

(9)

Calculate the difference between β₀ and β and take the absolute value:

$${\epsilon }_r=\left|\beta -{\beta }_0\right|$$

(55)

(10)

If ${\epsilon }_r>\epsilon$, where ε is the permissible error (e.g., $\epsilon ={10}^{-6}$), the new checking point in the original space can be obtained based on Equation (56):

$$x_1=x_0+J(u_1-u_0)$$

(56)

(11)

Substituting x₁ for x₀ in Step (3), repeat Steps (4) through (11) until convergence.

5. Discussion of the Applicable Bound of Equivalent Correlation Coefficients

For the novel structural reliability calculation method proposed in this paper, its applicability and scope of use need to be clarified. In this section, through mathematical derivation, detailed boundary formulas and range proofs are provided.

To simplify the description, Equation (36) was rewritten as follows:

$${\rho }_{ij}=A{\rho }_{0ij}^2+B{\rho }_{0ij}\begin{array}{cc}& \left(B>0\right)\end{array}$$

(57)

where

$$A=\frac{2c_i{c}_j}{{\sigma }_{Xi}{\sigma }_{Xj}}, B=\frac{b_i{b}_j}{{\sigma }_{Xi}{\sigma }_{Xj}}$$

(58)

According to Equations (25)–(27), the values of b_i (or b_j) will always be bigger than 0, resulting in a constant B > 0. The sign of A, however, is determined by the sign of c_ic_j, that is, λ_3Xi λ_3Xj.

If A ≠ 0, then Equation (57) can be expressed as a quadratic equation in terms of c_0ij, and is equivalent to:

$$\begin{array}{cc}\frac{{\rho }_{ij}}{A}={\rho }_{0ij}^2+\frac{B}{A}{\rho }_{0ij}& \left(B>0\right)\end{array}$$

(59)

To simplify the expression, the right-hand side of Equation (59) is expressed as h(ρ_0ij), i.e.,

$$\frac{{\rho }_{ij}}{A}=h\left({\rho }_{0ij}\right)={\rho }_{0ij}^2+\frac{B}{A}{\rho }_{0ij}$$

(60)

Subsequently, the extreme value of h(ρ_0ij) can be obtained as follows:

$$h_{ext}=h\left(\frac{-B}{2A}\right)=\frac{-B^2}{4A^2}$$

(61)

For A > 0, Equation (56) yields two real roots, namely ${\rho }_{0ij-1}=0$ and ${\rho }_{0ij-2}=-B/A<0$. Figure 1 depicts the shape of h(ρ_0ij) for A > 0, where the solid line denotes the region satisfying the condition that ${\rho }_{ij}\cdot {\rho }_{0ij}\ge 0$ and ρ_0ij is an increasing function of ρ_ij.

Based on Figure 1, the equivalent correlation coefficient ρ_0ij can be determined for A > 0 as follows:

$${\rho }_{0ij}=\frac{\sqrt{B^2+4A{\rho }_{ij}}-B}{2A}$$

(62)

To ensure that ρ_0ij obtained from Equation (62) satisfies the definition of correlation coefficient (i.e., $-1\le {\rho }_{0ij}\le 1$), an applicable bound for ρ_ij exists, i.e., ${\rho }_{ij}\in [{\rho }_{ij-\min },{\rho }_{ij-\max }]$, where ρ_ij-max and ρ_ij-min represent the upper and lower bounds, respectively.

When h_ext > −1/A (B² − 4A < 0), i.e., that is, the interval $\left(-\infty ,0\right]$ of ρ_0ij can be partitioned into two subintervals using the abscissas corresponding to ρ_0ij = −B/2A, as illustrated in Figure 1b. If ρ_0ij = −1 is in Interval I, the lower bound is determined as ${\rho }_{ij-\min }=-B^2/4A$. On the other hand, if ρ_0ij = −1 is in Interval II, the lower bound is ${\rho }_{ij-\min }=A-B$ to ensure ${\rho }_{0ij}\ge -1$. The lower bound for B²−4A < 0 can be expressed as ${\rho }_{ij-\min }=\max \left\{-B^2/4A,A-B\right\}$.

When h_ext ≤ −1/A (B² ≥ 4A), the interval $\left(-\infty ,0\right]$ of ρ_0ij can be partitioned into two subintervals using the abscissas corresponding to ρ_ij = −1, as depicted in Figure 1b. If ρ_0ij = −1 is in Interval I, then the lower bound is determined ${\rho }_{ij-\min }=A-B$ to ensure ${\rho }_{0ij}\ge -1$. Conversely, if ρ_0ij = −1 is in Interval II, the lower bound is set to ρ_ij-min = −1. For B² ≥ 4A, the lower bound can be expressed as ${\rho }_{ij-\min }=\max \left\{-1,A-B\right\}$.

Based on the previous discussion, if A > 0, the appropriate maximum value of the original correlation coefficient, which is equivalent to the lower bound of ρ_ij, can be summarized as follows:

$${\rho }_{ij-\min }=\{\begin{array}{l}\max \left\{-1,A-B\right\}, B^2-4A\ge 0\\ \max \left\{-B^2/4A,A-B\right\}, B^2-4A<0\end{array}$$

(63)

Substituting Equation (58) into the above Equation (63), the lower bound of ρ_ij can be expressed as follows:

$${\rho }_{ij-\min }=\{\begin{array}{l}\max \left\{-1,\left(2c_i{c}_j-b_i{b}_j\right)/{\sigma }_i{\sigma }_j\right\}, c_i,c_j\in \mathsf{\Omega }\\ \max \left\{\frac{-b_i^2{b}_j^2}{{\sigma }_i{\sigma }_j{c}_i{c}_j},\left(2c_i{c}_j-b_i{b}_j\right)/{\sigma }_i{\sigma }_j\right\}, \mathrm{otherwise}\end{array}$$

(64)

where

$$\mathsf{\Omega }=\left\{c_i,c_j|c_i{c}_j\le b_i^2{b}_j^2/8{\sigma }_i{\sigma }_j\right\}$$

(65)

(1)

When A > 0, the upper bound of ρ_ij can also be expressed as follows:

$${\rho }_{ij-\max }=\min \left\{1,A+B\right\}$$

(66)

then,

$${\rho }_{ij-\max }=\min \left\{1,\left(2c_i{c}_j+b_i{b}_j\right)/{\sigma }_i{\sigma }_j\right\}$$

(67)

(2)

When A< 0, ρ_0ij is determined by Equation (62), the applicable bound of ρ_ij can be expressed as:

$${\rho }_{ij-\max }=\{\begin{array}{l}\min \left\{1,A+B\right\}, B^2+4A\ge 0\\ \min \left\{-B^2/4A,A+B\right\}, \mathrm{otherwise}\end{array}$$

(68)

$${\rho }_{ij-\min }=\max \left\{-1,A-B\right\}$$

(69)

(3)

If A = 0, Equation (59) degenerates as a linear equation about ρ_0ij as follows:

$${\rho }_{0ij}=\frac{{\rho }_{ij}}{B}$$

(70)

The applicable bound of ρ_ij is given as follows:

$${\rho }_{ij-\min }=\max \left\{-1,-B\right\}$$

(71)

$${\rho }_{ij-\max }=\min \left\{1,B\right\}$$

(72)

6. Numerical Examples

6.1. Example 1: Computational Procedure with a Simplified Performance Function

Examine a cantilevered beam, as depicted in Figure 2, which is subjected to a uniformly distributed load (q) and possesses a beam length (l) of 4.0 m. The beam materials are characterized by Young’s modulus (E), while the inertia moment of the beam section is denoted by (I). The permissible deflection at the unattached extremity is 1/200 of the total beam length. As a result, the performance function can be formulated as follows:

$$G\left(\mathbf{X}\right)=\frac{l}{200}-\frac{q{l}^3}{EI}=0.02-8\frac{q}{EI}$$

(73)

The performance function, G(X), is articulated as a function encompassing three random variables: q, E, and I. The respective probability distributions and correlation matrix for these variables are presented in

Table 3. Random Variable Distributions and Correlation Matrix for Example 1.

Variable	Distribution	Mean Value	Coefficient of Variation	Correlation Matrix
q (N/m)	Lognormal	1000	0.2	$\left[\begin{array}{ccc}1.0& 0.2& 0.2\\ 0.2& 1.0& 0.2\\ 0.2& 0.2& 1.0\end{array}\right]$
E (N/m2)	Lognormal	2 × 1010	0.05
I (m4)	Lognormal	3.9025 × 10−5	0.1

.

6.1.1. Case 1: Structural Reliability with Correlated Variables

According to Equations (3)–(5) or Table 1, the first three L-moments of three random variables (q, E, and I) can be respectively given as follows:

$$\begin{gathered} {\lambda }_{1q}=1000, {\lambda }_{2q}=111.369, {\lambda }_{3q}=10.753 \\ {\lambda }_{1E}=2 \times {10}^{10}, {\lambda }_{2E}=5.637 \times {10}^8, {\lambda }_{3E} =1.376 \times {10}^7 \\ {\lambda }_{1I}=3.9025 \times {10}^{-5}, {\lambda }_{2I}=2.194\times {10}^{-6}, {\lambda }_{3I}=1.069\times {10}^{-7} \end{gathered}$$

Using Equations (25)–(27), the corresponding polynomial coefficients are given as follows:

$$\begin{gathered} a_1=980.389, b_1=197.396, c_1=19.504 \\ {a}_2=2 \times {10}^{10}, b_2=9.992 \times {10}^8, c_2=2.496 \times {10}^7 \\ {a}_3=3.883 \times {10}^{-5}, b_3=3.889 \times {10}^{-6}, c_3=1.937 \times {10}^{-7} \end{gathered}$$

Using Equation (58), the parameter matrix A can be obtained as follows:

$$A=\left[\begin{array}{ccc}0.0190199& 0.00486773& 0.00969043\\ 0.00486773& 0.00124579& 0.00248006\\ 0.00969043& 0.00248006& 0.00493717\end{array}\right]$$

Upon observation, it can be seen that all elements of A are greater than 0. Subsequently, by utilizing Equations (64) and (67), the upper and lower bounds of ρ_ij can be expressed as follows:

$$\begin{gathered} {\rho }_{ij-\max }=\left[\begin{array}{ccc}0.993099& 0.990987& 0.993189\\ 0.990913& 0.999532& 0.998103\\ 0.993198& 0.998095& 0.997911\end{array}\right] \\ {\rho }_{ij-\min }=\left[\begin{array}{ccc}-0.955019& -0.981263& -0.973075\\ -0.981249& -0.997102& -0.993372\\ -0.973798& -0.992997& -0.987993\end{array}\right] \end{gathered}$$

The pertinent range of initial correlation coefficients matrix can be expeditiously procured as $\Theta =\left[-0.955019, 0.990913\right]$, which infers that the initial correlation matrix presented in this instance is suitable for the third L-moments transformation method. It is evident that the applicable range of the original correlation coefficient is sufficiently expansive, thereby rendering it satisfactory for a diverse array of general engineering problems.

Based on Equation (62), the equivalent correlation matrix is given as follows:

$${\mathbf{C}}_{\mathbf{Z}}=\left[\begin{array}{ccc}1.0000& 0.2019& 0.2020\\ 0.2019& 1.0000& 0.2005\\ 0.2020& 0.2005& 1.0000\end{array}\right]$$

The Cholesky decomposition results of C_Z can then be expressed as follows:

$$L_0=\left[\begin{array}{ccc}1.0& 0& 0\\ 0.2019& 0.9794& 0\\ 0.2020& 0.1631& 0.9657\end{array}\right], L_0^{-1}=\left[\begin{array}{ccc}1.0& 0& 0\\ -0.2061& 1.0210& 0\\ -0.1744& -0.1724& 1.0355\end{array}\right]$$

Assume initial checking point as:

$$x_0={\left(q, E, I\right)}^{\mathrm{T}}={\left(1000, 2\times {10}^{10}, 3.9025\times {10}^{-5}\right)}^{\mathrm{T}}$$

The corresponding checking point in standard normal space can be readily obtained as $u_0={\left(0.0981,0.0052,0.0298\right)}^{\mathrm{T}}$ using Equation (50). Subsequently, the initial reliability index β₀ at u₀ can be calculated as ${\beta }_0=\sqrt{u_0^{\mathrm{T}}{u}_0}=0.1025$.

The Jacobian matrix at u₀ can be evaluated using Equation (51) as follows:

$$J=\left(\begin{array}{ccc}201.214& 0& 0\\ 2.0193\times {10}^8& 9.7979\times {10}^8& 0\\ 7.8954\times {10}^{-7}& 6.3739\times {10}^{-7}& 3.7740\times {10}^{-6}\end{array}\right)$$

Using Equation (52), the value of the performance function and gradient vector at the checking point u₀ can be evaluated as follows:

$$\begin{gathered} g\left(u_0\right)=G\left(x_0\right)=0.00969 \\ \nabla g(u_0)=\left(\begin{array}{c}-1.752\times {10}^{-3}\\ 6.695\times {10}^{-4}\\ 9.912\times {10}^{-4}\end{array}\right) \end{gathered}$$

Then, the next checking point u₁ can be obtained based on Equation (53) as:

$$u_1=\frac{1}{{\nabla }^{\mathrm{T}}g(u_0)\cdot \nabla g(u_0)}\nabla g(u_0)\left[{\nabla }^{\mathrm{T}}g(u_0)\cdot u_0-g(u_0)\right]=\left(\begin{array}{c}3.8499\\ -1.4713\\ -2.1791\end{array}\right)$$

Then, the reliability index can be obtained as $\beta =\sqrt{u_1^T{u}_1}=4.66195$.

The absolute difference between β₀ and β can be determined as:

$${\epsilon }_r=\frac{|\beta -{\beta }_0|}{\beta }=0.9780>1.0\times {10}^{-6}$$

Checking point in original space is determined as:

$$x_1=x_0+J(u_1-u_0)={(1754.95,1.93105\times {10}^{10},0.0000327)}^T$$

The above procedure is repeated using x₁ as the new checking point until convergence is achieved. A comparison of results from different methods for Case 1 is presented in

Table 4. Comparison of results by different methods for Case 1.

Method	Design Point u*	Design Point x*	β	Pf
L-Moments	(2.7039, −1.1539, −1.7016)T	(1627.03, 1.9415 × 1010, 3.3415 × 10−5)T	3.395	3.40 × 10−4
C-Moments	(2.7031, −1.1527, −1.6949)T	(1633.32, 1.9376 × 1010, 3.3763 × 10−5)T	3.391	3.47 × 10−4
MCS	--	--	3.397	3.36 × 10−4

. The convergence criterion is met in five steps, resulting in a reliability index of β_LM = 3.395 at the design point u* = (2.7039, −1.1539, −1.7016)^T in standard normal space, which corresponds to the design point x* = (1627.03, 1.9415 × 10¹⁰, 3.3415 × 10⁻⁵)^T in the original space.

A comparative reliability analysis was conducted using the first three C-moments to perform the x-u and u-x transformations for the same example. The resulting reliability index, denoted as β_CM, was obtained as 3.391 after convergence was achieved (${\epsilon }_r\le 1.0\times {10}^{-6}$) in five steps. The corresponding design point in standard normal space was u* = (2.7031, −1.1527, −1.6949)^T, and the corresponding design point in the original space was x* = (1633.32, 1.9376 × 10¹⁰, 3.3763 × 10⁻⁵)^T.

Based on the available distribution information and initial correlation matrix of the random variables, the failure probability is determined to be 3.36 × 10⁻⁴ with a coefficient of variation of 3.77%, using the Monte Carlo simulation (MCS) approach with a sample size of 2 × 10⁶. The corresponding reliability index (β_MCS) is found to be 3.397. Notably, the results obtained from the first three L-moments- and C-moments-based FORM approaches are found to be in close agreement with those obtained from MCS. It is noteworthy that the reliability index and failure probability ascertained via the L-moments method exhibit greater proximity to the values computed using the MCS method. This observation suggests that the L-moments approach yields superior accuracy in comparison to the C-moments method.

6.1.2. Case 2: Examining Random Variables with Greater Skewness

The proposed method, premised on the first three L-moments in Equation (16) and Equations (25)–(27), circumvents limitations on L-skewness, in contrast to the conventional third-moment standardization. If the coefficient of variation of q is 0.8, its skewness attains 2.924, which surpasses the scope of the third-moment transformation predicated on C-moments. In this scenario, the initial correlation matrix is verified with $\Theta =\left[-0.539, 0.861\right]$, facilitating reliability analysis based on L-moments. The acquired reliability index β_LM = 1.266, converges to u* = (1.2429, −0.1334, −0.1972)^T, corresponding to x* = (1963.9396, 2.0201 × 10¹⁰, 3.9149 × 10⁻⁵)^T, after reaching convergence (${\epsilon }_r\le 1.0\times {10}^{-6}$) in five steps. Utilizing MCS, the failure probability amounts to 9.272 × 10⁻² (the coefficient of variation of MCS is 0.991%), with the corresponding reliability index of 1.319. It can be observed that FORM, based on L-moments, yields a result comparable to that of MCS, while FORM, based on C-moments, cannot be executed for this case.

Figure 3 illustrates the reliability indices in relation to the coefficients of variation of q, as determined using FORM based on L-moments and C-moments, along with the MCS results. Upon examining Figure 3, it becomes evident that the L-moments-based FORM results consistently demonstrate superior accuracy and alignment with the MCS outcomes, as compared to the C-moments-based FORM results, when the coefficients of variation of q range from 0.2 to 0.778.

As the coefficients of variation of q exceed 0.778, the L-moments-based FORM continues to provide comparable results to those obtained through MCS, further emphasizing its advantageous performance. In contrast, the C-moments-based FORM becomes impractical, as the skewness of q surpasses the application range of the third-moment transformation predicated on C-moments. This observation accentuates the superior performance of the L-moments method compared to the C-moments method in the given context.

6.2. Example 2: Reliability Analysis with Unknown Probability Distribution

The second example encompasses a roof structure subjected to a uniformly distributed vertical load (q_r), as delineated in Figure 4, which has been previously examined by [2]. The top cords and compression bars are composed of concrete, while the bottom cords and tension bars are constituted of steel. In the structural analysis, the uniformly distributed load (q_r) was transformed into three nodal loads, each denoted as P = q_rl/4. The serviceability limit state of the structure, with respect to its maximum vertical displacement, was considered. The performance function is expressed as follows:

$$G\left(X\right)=0.035-\frac{q{l}^2}{2}(\frac{3.81}{A_c{E}_c}+\frac{1.13}{A_s{E}_s})$$

(74)

where l represents the roof span, and the admissible deflection is arbitrarily fixed to 0.035 m. Furthermore, (E_s, A_s) and (E_c, A_c) correspond to the Young’s modulus and cross-sectional areas of the steel and concrete beams, respectively.

Table 5. Random Variable Distributions and Sample L-moments for Example 2.

Variables	Distributions	Mean Value	C.O.V	First Three Sample L-Moments
				l1	l2	l3
qr (N/m)	Gumbel	2 × 104	0.10	2000.85	1081.72	0.1707
l (m)	Gumbel	12	0.01	11.99	0.0648	0.1698
As (m2)	Normal	9.82 × 10−4	0.06	9.83 × 10−4	3.31 × 10−5	3.35 × 10−4
Ac (m2)	Normal	400 × 10−4	0.12	3.99 × 10−4	2.71 × 10−3	3.15 × 10−4
Es (N/m2)	Lognormal	1 × 1011	0.06	9.99 × 1010	3.37 × 109	0.0292
Ec (N/m2)	Lognormal	2 × 1010	0.06	1.99 × 1010	6.77 × 109	0.0297

consolidates the statistical information pertaining to the random variables.

Assuming all random variables possess known CDFs/PDFs, the reliability index and failure probability can be effectively ascertained through the FORM procedure and MCS method. In order to examine the efficiency of the proposed approach, which incorporates random variables with unknown CDFs/PDFs, the CDFs/PDFs of random variables q_r, l, A_s, A_c, E_s, and E_c are considered unknown, with only a limited number of measured samples assumed to be known. Consequently, the first three sample L-moments can be calculated employing Equations (12)–(14) and Equation (15) and are presented in the rightmost section of Table 5. By utilizing the first three L-moments, the u-x and x-u transformations of the random variables can be seamlessly executed in accordance with Section 3. Subsequently, the reliability index and failure probability including random variables with unknown CDFs/PDFs, can also be readily obtained.

Upon examination of the available distribution information for the random variables, the failure probability is ascertained to be 1.96 × 10⁻⁴, accompanied by a coefficient of variation equal to 2.87%. This is achieved through the application of the MCS method, employing a sample size of 5 × 10⁶. The corresponding reliability index (β_MCS) is determined to be 3.542. With respect to the known distribution information of random variables, the first-order reliability index (β_FORM) is calculated as 3.679 using a general FORM procedure. In instances where distribution information remains unknown, the reliability index (β_LM) computed by amalgamating FORM with merely the first three sample L-moments is 3.594.

It is worth noting that the results procured from the FORM approaches only based on the first three sample L-moments exhibit close alignment with those derived from MCS and the general FORM procedure. The comprehensive results acquired through the utilization of the CDFs/PDFs and the first three sample L-moments of q_r, l, A_s, A_c, E_s, and E_c are presented in

Table 6. Comparisons of FORM Procedure Using Different Methods for Example 2.

Variables	Using CDFs/PDFs			Using Sample L-Moments
	Checking Point		Jacobian	Checking Point		Jacobian
	x*	u*	du/dx	x*	u*	du/dx
qr (N/m)	22423.79	1.343	1379.57	22447.15	1.334	1381.84
l (m)	12.09	0.432	0.119	11.97	0.479	0.122
As (m2)	9.24 × 10−4	−0.981	5.78 × 10−5	9.23 × 10−4	−0.979	5.82 × 10−5
Ac (m2)	0.037	−0.863	0.0049	0.036	−0.861	0.0047
Es (N/m2)	9.47 × 1010	−0.991	5.97 × 109	9.33 × 1010	−0.997	5.97 × 109
Ec (N/m2)	1.94 × 1010	−0.402	1.19 × 109	1.94 × 1010	−0.405	1.19 × 109

. The table demonstrates that the checking points (in both original and standard normal spaces), the Jacobian for the last iteration, employing the first three sample L-moments of q_r, l, A_s, A_c, E_s, and E_c, exhibit a general proximity to those ascertained using the CDFs/PDFs of q_r, l, A_s, A_c, E_s, and E_c. This observation indicates that the first-order reliability index and calculation procedure, derived using exclusively the first three sample L-moments of random variables, exhibit substantial coherence with the index and calculation process ascertained via the CDFs/PDFs of random variables. This evidence not only reinforces the feasibility of the method proposed in this paper when the distribution of random variables is unknown but also accentuates its alignment with the scenarios encountered in engineering practice. Consequently, it encourages the promotion and application of the proposed method within real-world engineering contexts.

6.3. Example 3: Reliability Analysis for Implicit Performance Functions

Consider a two-bay five-story frame structure subjected to lateral loads as depicted in Figure 5. The cross-sectional areas and moments of inertia of the frame members are listed in

Table 7. Probabilistic Characteristics of the Basic Random Variables.

Variable	Distribution	The First Three L-Moments			The Corresponding Parameters
		λ1	λ2	λ3	a	b	c
P1	Weibull	135	15.213	−0.9698	136.759	26.9643	−1.75902
P2	Weibull	90	10.142	−0.6465	91.1726	17.9762	−1.17262
P3	Weibull	70	7.888	−0.5029	70.9122	13.9811	−0.91216
E1	Lognormal	2.1 × 107	591906	14449.3	2.10 × 107	1.05 × 106	26208.1
E2	Lognormal	2.4 × 107	676464	16513.5	2.40 × 107	1.20 × 106	29952.2
I1	Lognormal	8.1 × 10−3	4.6 × 10−4	2.2 × 10−5	8.07 × 10−3	8.15 × 10−4	3.99 × 10−5
I2	Lognormal	0.011	6.2 × 10−4	3.0 × 10−5	0.0109	0.0011	5.44 × 10−5
I3	Lognormal	0.0213	0.0012	5.8 × 10−5	0.0212	0.00213	1.05 × 10−4
I4	Lognormal	0.0295	0.0015	7.1 × 10−5	0.0258	0.00266	0.00013
I5	Lognormal	0.0108	6.1 × 10−4	2.9 × 10−5	0.0107	0.00108	5.26 × 10−5
I6	Lognormal	0.0141	7.9 × 10−4	3.9 × 10−5	0.0140	0.00140	7.07 × 10−5
I7	Lognormal	0.0232	1.3 × 10−3	6.4 × 10−5	0.0231	0.00230	1.16 × 10−4
A1	Lognormal	0.312	0.0175	8.5 × 10−4	0.3105	0.03102	0.00154
A2	Lognormal	0.372	0.0209	1.0 × 10−3	0.3702	0.03704	0.00181
A3	Lognormal	0.505	0.0284	1.4 × 10−3	0.5025	0.05034	0.00254
A4	Lognormal	0.557	0.0313	1.5 × 10−3	0.5543	0.05548	0.00272
A5	Lognormal	0.253	0.0142	6.9 × 10−4	0.2517	0.02517	0.00125
A6	Lognormal	0.291	0.0164	8.0 × 10−4	0.2895	0.02907	0.00145
A7	Lognormal	0.372	0.0209	1.0 × 10−3	0.3702	0.0370	0.00181

. The performance function is given as follows:

$$G(X)=u_{\lim }-u(X)$$

(75)

The specified limit u_lim = 0.055 m and the top-floor displacement u(X), which is determined by finite element analysis, are used in the following performance function. The correlation matrix of basic random variables is presented below:

$$C_X=\begin{array}{c}{P}_1\\ P_2\\ P_3\\ E_1\\ E_2\\ I_1\\ ⋮\\ A_7\end{array}\left(\begin{array}{cccccccc}1& & & & & & & \\ 0.5& 1& & & & & & \\ 0.5& 0.5& 1& & & & & \\ 0& 0& 0& 1& & & & \\ 0& 0& 0& 0.9& 1& & & \\ 0& 0& 0& 0& 0.13& 1& & \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & \ddots & \\ 0& 0& 0& 0& 0.13& \cdots & 0.13& 1\end{array}\right)$$

(76)

The probabilistic characteristics and associated parameters of the random variables are presented in Table 7. The transformations from the standard normal space to the original space and vice versa are conducted for the random variables using Equations (45) and (50), respectively.

In this case, the boundary of the initial correlation matrix is examined to satisfy $\Theta =\left[-0.981,0.985\right]$, thereby enabling the implementation of reliability analysis based on L-moments.

For the implicit performance function in this case, numerical differentiation methods can be utilized to determine the gradient vector $\nabla G(x_0)$. If the central difference method is employed, the expression for $\nabla G(x_0)$ is given by:

$$\nabla G(x_0)=\left(\begin{array}{c}\frac{G(x_{01}+h,x_{02},\cdots ,x_{0n})-G(x_{01}-h,x_{02},\cdots ,x_{0n})}{2h}\\ ⋮\\ \frac{G(x_{01},\cdots ,x_{0i}+\Delta ,\cdots ,x_{0n})-G(x_{01},\cdots ,x_{0i}-\Delta ,\cdots ,x_{0n})}{2\Delta }\\ ⋮\\ \frac{G(x_{01},x_{02},\cdots ,x_{0n}+\Delta )-G(x_{01},x_{02},\cdots ,x_{0n}-\Delta )}{2\Delta }\end{array}\right)$$

(77)

where x_0i (i = 1,2, …, n) represents the ith element of vector x₀; and $\Delta ={10}^{-4}$ is the step size (generally chosen to be small enough as required).

Once the convergence is achieved (with ${\epsilon }_r\le 1.0\times {10}^{-4}$) in five steps, the reliability index can be easily determined as β_LM = 2.209. Moreover, employing Monte Carlo simulation with known PDFs/CDFs and the initial correlation matrix, the reliability index can be calculated as β_MCS = 2.215 (with corresponding failure probability P_f = 1.337 × 10⁻²) with high accuracy (coefficient of variation for P_f = 2.721%). Thus, the results obtained from the two methods in this example are almost identical.

7. Concluding Remarks

This study proposes a second-order polynomial normal transformation for structural reliability analysis based on the first three L-moments for both independent and correlated variables. The criteria for the admissible range of the initial correlation matrix have also been explored. Moreover, a FORM method has been developed based on the first three L-moments. Numerical examples are presented to investigate the accuracy and efficiency of the proposed method for structural reliability analysis. Several conclusions are drawn, as summarized in the following:

(1)

In practical engineering, the probability information of random variables is often incomplete, with access to only limited statistical data, rendering reliability calculations unfeasible. To address this issue, this paper presents an innovative solution. The proposed solution consists of two steps: first, employing the normal transformation and inverse transformation formulas based on the first three L-moments proposed in this paper, random variables can be approximated using polynomials containing standard normal random variables. Subsequently, after obtaining the approximate expressions for the random variables, the traditional FORM procedure is combined to calculate the reliability of the performance function using the approximate distribution and correlation matrix of the random variables. Several numerical examples provided, compared with the MCS method, demonstrate the feasibility and accuracy of the proposed method.

(2)

Although the normal transformation and inverse transformation of random variables can be achieved using C-moments, the accuracy of the C-moments-based FORM method tends to degrade as the coefficient of variation increases. When the coefficient of variation exceeds 0.778, the C-moments method becomes completely ineffective. In contrast, fortunately, the L-moments-based FORM method proposed in this paper consistently maintains agreement with the reliability indices calculated using the MCS method, regardless of the magnitude of the coefficient of variation. This is vividly illustrated in Case 2 of Example 2, which further demonstrates the broader applicability of the proposed method without being constrained by the coefficient of variation. It is worth mentioning that this paper also provides detailed boundary formulas and range justification for the equivalent correlation coefficients of the proposed method through mathematical derivation, further clarifying the applicability and scope of the new method.

(3)

In this study, a total of three numerical examples have been investigated to substantiate the feasibility, accuracy, and simplicity of the method proposed herein. These examples are grounded in real-world scenarios, encompassing a wealth of engineering contexts. They have been carefully selected to address diverse distribution types, multiple random variables, and intricate structural configurations from a variety of perspectives, underlining the extensive applicability of the presented method. Furthermore, the derived normal transformation and inverse transformation formulas for random variables, based on L-moments, are notably concise and straightforward. This simplicity not only facilitates comprehension among engineering practitioners but also promotes ease of implementation in real-life engineering problems. In light of these findings, it is recommended that the method proposed in this paper be further disseminated and applied within the engineering community. By doing so, it is anticipated that the method will contribute significantly to the advancement of structural reliability analysis and offer valuable insights for practitioners seeking to address complex engineering challenges.

Regarding future work, two primary research directions can be pursued to build upon the findings of this study. Firstly, a deeper investigation of the proposed method’s applicability to dynamic reliability analysis and more complex structural systems is needed. By extending the method to these domains, its versatility and robustness can be further established. Secondly, the exploration of alternative techniques that can substitute for normal transformation and inverse transformation is of interest, such as deep learning technologies. The potential integration of advanced machine learning algorithms with the proposed method could lead to enhanced performance and new insights into the field of structural reliability analysis.

By addressing these research avenues, the proposed method can continue to evolve and contribute to the advancement of structural reliability analysis and the resolution of complex engineering challenges.