Second order reliability method

For the second order reliability method, the random variables \(X\) first should be transformed to standard normal uncorrelated variables \(U\). Then, the limit state function \(g(U)\) can be approximated with a second order Taylor expansion,

\[g(U) \approx g(U^*) + \nabla g(U^*)^T(U-U^*) + \frac{1}{2}(U-U^*)^T \nabla^2 g(U^*) (U-U^*)\]

where \(U^*\) is the design point or most probable failure point (MPP); \(\nabla^2 g(U^*)\) represent the Hessian matrix evaluated at the design point, it can be represented by,

\[\nabla^2 g(U^*)_{ij} = \frac{\partial g(U^*)}{\partial u_i \partial u_j}\]

Since \(g(U^*)=0\), the second order Taylor expansion can be represented by,

\[g(U) \approx \nabla g(U^*)^T(U-U^*) + \frac{1}{2}(U-U^*)^T \nabla^2 g(U^*) (U-U^*)\]

To solve the problem, a transformation \(Y=HU\) is performed so that the last coordinate coincides with the vector \(U^*\) from the origin to the design point (\(\beta\) vector). \(H\) can be obtained by a Gram-Schmidt orthogonalization. Thus, the Taylor expansion is,

\[g(Y) \approx -y_n + \beta + \frac{1}{2}(Y-Y^*)^T H \frac{\nabla^2 g(U^*)}{||\nabla g(U^*) ||} H^T (Y-Y^*)\]

where \(Y^* = \{0,0,\dots, \beta\}^T\) is the design point in \(Y\) space corresponding to the design point \(U^*\) in \(U\) space; \((Y-Y^*) = \{y_1,y_2,\dots, y_n-\beta\}^T\).

The main curvatures \(k_i\) can be obtained by transforming the \((n-1)\times(n-1)\) order matrix of \(H \frac{\nabla^2 g(U^*)}{||\nabla g(U^*) ||} H^T\) to a diagonal matrix (i.e., eigenvalues).

References:

• Choi, S.K., Canfield, R.A. and Grandhi, R.V., 2007. Reliability-Based Structural Design. Springer London.

# Import auxiliary libraries for demonstration

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

plt.rcParams[ "figure.figsize" ] = [ 5, 4 ]

plt.rcParams[ "figure.dpi" ] = 80
plt.rcParams[ "font.family" ] = "Times New Roman"
plt.rcParams[ "font.size" ] = '14'

Breitung SORM

With the \(n-1\) main curvatures \(k_i\) and first order reliability index \(\beta\), Breitung proposed the parabolic approximation of probability of failure,

\[p_f = \Phi(-\beta) \prod_{i=1}^{n-1} (1+k_i \beta)^{-1/2}\]

where \(\Phi(\cdot)\) is the standard normal Cumulative Distribution Function (CDF); \(\beta\) is the first order reliability index; \(k_i\) is the main curvatures of the limit-state function at design point; It should be noted that Breitung formula is applicable for large \(\beta\).

Function breitungSORM implements the SORM with Breitung algorithm. The Nataf transformation is used in the method to map the random variables from X space to U space.

References:

• Breitung, K., 1984. Asymptotic approximations for multinormal integrals. Journal of Engineering Mechanics, 110(3), pp.357-366.

• Hu, Z., Mansour, R., Olsson, M. and Du, X., 2021. Second-order reliability methods: a review and comparative study. Structural and Multidisciplinary Optimization, 64(6), pp.3233-3263.

• Bourinet, J.M., 2018. Reliability analysis and optimal design under uncertainty-Focus on adaptive surrogate-based approaches (Doctoral dissertation, Université Clermont Auvergne).

Function help

from ffpack.rrm import breitungSORM
help( breitungSORM )

Help on function breitungSORM in module ffpack.rrm.secondOrderReliabilityMethod:

breitungSORM(dim, g, dg, distObjs, corrMat, quadDeg=99, quadRange=8, dx=1e-06)
    Second order reliability method based on Breitung algorithm.

    Parameters
    ----------
    dim: integer
        Space dimension ( number of random variables ).
    g: function
        Limit state function. It will be called like g( [ x1, x2, ... ] ).
    dg: array_like of function
        Gradient of the limit state function. It should be an array_like of function
        like dg = [ dg_dx1, dg_dx2, ... ]. To get the derivative of i-th random
        variable at ( x1*, x2*, ... ), dg[ i ]( x1*, x2*, ... ) will be called.
        dg can be None, see the following Notes.
    distObjs: array_like of distributions
        Marginal distribution objects. It should be the freezed distribution
        objects with pdf, cdf, ppf. We recommend to use scipy.stats functions.
    corrMat: 2d matrix
        Correlation matrix of the marginal distributions.
    quadDeg: integer
        Quadrature degree for Nataf transformation
    quadRange: scalar
        Quadrature range for Nataf transformation. The integral will be performed
        in the range [ -quadRange, quadRange ].
    dx : scalar, optional
        Spacing for auto differentiation. Not required if dg is provided.

    Returns
    -------
    beta: scalar
        Reliability index.
    pf: scalar
        Probability of failure.
    uCoord: 1d array
        Design point coordinate in U space.
    xCoord: 1d array
        Design point coordinate in X space.

    Raises
    ------
    ValueError
        If the dim is less than 1.
        If the dim does not match the disObjs and corrMat.
        If corrMat is not 2d matrix.
        If corrMat is not positive definite.
        If corrMat is not symmetric.
        If corrMat diagonal is not 1.

    Notes
    -----
    If dg is None, the numerical differentiation will be used. The tolerance of the
    numerical differentiation can be changed in globalConfig.

    Examples
    --------
    >>> from ffpack.rrm import breitungSORM
    >>> dim = 2
    >>> g = lambda X: -np.sum( X ) + 1
    >>> dg = [ lambda X: -1, lambda X: -1 ]
    >>> distObjs = [ stats.norm(), stats.norm() ]
    >>> corrMat = np.eye( dim )
    >>> beta, pf, uCoord, xCoord = breitungSORM( dim, g, dg, distObjs, corrMat )

Example with explicit derivative of LSF

# Define the dimension for the FORM problem
dim = 2

# Define the limit state function (LSF) g
g = lambda X: X[ 0 ] ** 4 + 2 * X[ 1 ] ** 4 - 20

# Explicit derivative of LSF
# dg is a list in which each element is a partial derivative function of g w.r.t. X
# dg[0] = partial g / partial X[0]
# dg[1] = partial g / partial X[1]
dg = [ lambda X: 4 * X[ 0 ] ** 3, lambda X: 8 * X[ 1 ] ** 3 ]

# Marginal distributions and correlation Matrix of the random variables
distObjs = [ stats.norm( 5.0, 5.0 ), stats.norm( 5.0, 5.0 ) ]
corrMat = np.eye( dim )

beta, pf, uCoord, xCoord = breitungSORM( dim, g, dg, distObjs, corrMat )

print( "Reliability index: " )
print( beta )
print()
print( "Failure probability: " )
print( pf )
print()
print( "Design point coordinate in U space: " )
print( uCoord )
print()
print( "Design point coordinate in X space: " )
print( xCoord )

Reliability index:
0.9519628114174661

Failure probability:
0.06374492261292942

Design point coordinate in U space:
[-0.6398897170780924, -0.7048222075811494]

Design point coordinate in X space:
[1.8005514146095387, 1.4758889620942535]

Example with automatic differentiation of LSF

# Define the dimension for the FORM problem
dim = 2

# Define the limit state function (LSF) g
g = lambda X: X[ 0 ] ** 4 + 2 * X[ 1 ] ** 4 - 20

# If dg is None, the internal automatic differentiation algorithm will be used
dg = None

# Marginal distributions and correlation Matrix of the random variables
distObjs = [ stats.norm( 5.0, 5.0 ), stats.norm( 5.0, 5.0 ) ]
corrMat = np.eye( dim )

beta, pf, uCoord, xCoord = breitungSORM( dim, g, dg, distObjs, corrMat )

print( "Reliability index: " )
print( beta )
print()
print( "Failure probability: " )
print( pf )
print()
print( "Design point coordinate in U space: " )
print( uCoord )
print()
print( "Design point coordinate in X space: " )
print( xCoord )

Reliability index:
0.9519628114174661

Failure probability:
0.06374492260984004

Design point coordinate in U space:
[-0.6398897170780924, -0.7048222075811494]

Design point coordinate in X space:
[1.8005514146095387, 1.4758889620942535]

Tvedt SORM

Tvedt further derived a three-term approximation by ignoring terms of orders higher than two,

\[T_1 = \Phi(-\beta) \prod_{i=1}^{n-1} (1+k_i \beta)^{-1/2}\]

\[T_2 = \left[ \beta \Phi(-\beta) - \phi(\beta) \right] \left[ \prod_{i=1}^{n-1} (1+k_i \beta)^{-1/2} - \prod_{i=1}^{n-1} (1+k_i (\beta+1))^{-1/2} \right]\]

\[T_3 = (\beta + 1) \left[ \beta \Phi(-\beta) - \phi(\beta) \right] \left[ \prod_{i=1}^{n-1} (1+k_i \beta)^{-1/2} - \text{Re} \left[ \prod_{i=1}^{n-1} (1+k_i (\beta+1))^{-1/2} \right] \right]\]

\[p_f = T_1 + T_2 + T_3\]

where \(\Phi(\cdot)\) is the standard normal Cumulative Distribution Function (CDF); \(\phi(\cdot)\) is the standard normal Probabilistic Distribution Function (PDF); \(\beta\) is the first order reliability index; \(k_i\) is the main curvatures of the limit-state function at design point; \(\text{Re}\) is the real part of a complex number. It can be found that the \(T_1\) is the Breitung’s equation. The \(T_2\) and \(T_3\) can be interpreted as the correctors to the Breitung’s formula to increase the accuracy for moderate values of \(\beta\).

References:

• Tvedt, L., 1990. Distribution of quadratic forms in normal space—application to structural reliability. Journal of engineering mechanics, 116(6), pp.1183-1197.

• Hu, Z., Mansour, R., Olsson, M. and Du, X., 2021. Second-order reliability methods: a review and comparative study. Structural and Multidisciplinary Optimization, 64(6), pp.3233-3263.

• Bourinet, J.M., 2018. Reliability analysis and optimal design under uncertainty-Focus on adaptive surrogate-based approaches (Doctoral dissertation, Université Clermont Auvergne).

Function help

from ffpack.rrm import tvedtSORM
help( tvedtSORM )

Help on function tvedtSORM in module ffpack.rrm.secondOrderReliabilityMethod:

tvedtSORM(dim, g, dg, distObjs, corrMat, quadDeg=99, quadRange=8, dx=1e-06)
    Second order reliability method based on Tvedt algorithm.

    Parameters
    ----------
    dim: integer
        Space dimension ( number of random variables ).
    g: function
        Limit state function. It will be called like g( [ x1, x2, ... ] ).
    dg: array_like of function
        Gradient of the limit state function. It should be an array_like of function
        like dg = [ dg_dx1, dg_dx2, ... ]. To get the derivative of i-th random
        variable at ( x1*, x2*, ... ), dg[ i ]( x1*, x2*, ... ) will be called.
        dg can be None, see the following Notes.
    distObjs: array_like of distributions
        Marginal distribution objects. It should be the freezed distribution
        objects with pdf, cdf, ppf. We recommend to use scipy.stats functions.
    corrMat: 2d matrix
        Correlation matrix of the marginal distributions.
    quadDeg: integer
        Quadrature degree for Nataf transformation
    quadRange: scalar
        Quadrature range for Nataf transformation. The integral will be performed
        in the range [ -quadRange, quadRange ].
    dx : scalar, optional
        Spacing for auto differentiation. Not required if dg is provided.

    Returns
    -------
    beta: scalar
        Reliability index.
    pf: scalar
        Probability of failure.
    uCoord: 1d array
        Design point coordinate in U space.
    xCoord: 1d array
        Design point coordinate in X space.

    Raises
    ------
    ValueError
        If the dim is less than 1.
        If the dim does not match the disObjs and corrMat.
        If corrMat is not 2d matrix.
        If corrMat is not positive definite.
        If corrMat is not symmetric.
        If corrMat diagonal is not 1.

    Notes
    -----
    If dg is None, the numerical differentiation will be used. The tolerance of the
    numerical differentiation can be changed in globalConfig.

    Examples
    --------
    >>> from ffpack.rrm import tvedtSORM
    >>> dim = 2
    >>> g = lambda X: -np.sum( X ) + 1
    >>> dg = [ lambda X: -1, lambda X: -1 ]
    >>> distObjs = [ stats.norm(), stats.norm() ]
    >>> corrMat = np.eye( dim )
    >>> beta, pf, uCoord, xCoord = tvedtSORM( dim, g, dg, distObjs, corrMat )

Example with explicit derivative of LSF

# Define the dimension for the FORM problem
dim = 2

# Define the limit state function (LSF) g
g = lambda X: X[ 0 ] ** 4 + 2 * X[ 1 ] ** 4 - 20

# Explicit derivative of LSF
# dg is a list in which each element is a partial derivative function of g w.r.t. X
# dg[0] = partial g / partial X[0]
# dg[1] = partial g / partial X[1]
dg = [ lambda X: 4 * X[ 0 ] ** 3, lambda X: 8 * X[ 1 ] ** 3 ]

# Marginal distributions and correlation Matrix of the random variables
distObjs = [ stats.norm( 5.0, 5.0 ), stats.norm( 5.0, 5.0 ) ]
corrMat = np.eye( dim )

beta, pf, uCoord, xCoord = tvedtSORM( dim, g, dg, distObjs, corrMat )

print( "Reliability index: " )
print( beta )
print()
print( "Failure probability: " )
print( pf )
print()
print( "Design point coordinate in U space: " )
print( uCoord )
print()
print( "Design point coordinate in X space: " )
print( xCoord )

Reliability index:
0.9519628114174661

Failure probability:
0.03604492191838041

Design point coordinate in U space:
[-0.6398897170780924, -0.7048222075811494]

Design point coordinate in X space:
[1.8005514146095387, 1.4758889620942535]

Example with automatic differentiation of LSF

# Define the dimension for the FORM problem
dim = 2

# Define the limit state function (LSF) g
g = lambda X: X[ 0 ] ** 4 + 2 * X[ 1 ] ** 4 - 20

# If dg is None, the internal automatic differentiation algorithm will be used
dg = None

# Marginal distributions and correlation Matrix of the random variables
distObjs = [ stats.norm( 5.0, 5.0 ), stats.norm( 5.0, 5.0 ) ]
corrMat = np.eye( dim )

beta, pf, uCoord, xCoord = tvedtSORM( dim, g, dg, distObjs, corrMat )

print( "Reliability index: " )
print( beta )
print()
print( "Failure probability: " )
print( pf )
print()
print( "Design point coordinate in U space: " )
print( uCoord )
print()
print( "Design point coordinate in X space: " )
print( xCoord )

Reliability index:
0.9519628114174661

Failure probability:
0.03604492191636101

Design point coordinate in U space:
[-0.6398897170780924, -0.7048222075811494]

Design point coordinate in X space:
[1.8005514146095387, 1.4758889620942535]

Hohenbichler and Rackwitz SORM

Hohenbichler and Rackwitz derived a closed form expression based on a Taylor expansion,

\[p_f = \Phi(-\beta) \prod_{i=1}^{n-1} \left(1+k_i \frac{\phi(\beta)}{\Phi(\beta)} \right)^{-1/2}\]

where \(\Phi(\cdot)\) is the standard normal Cumulative Distribution Function (CDF); \(\phi(\cdot)\) is the standard normal Probabilistic Distribution Function (PDF); \(\beta\) is the first order reliability index; \(k_i\) is the main curvatures of the limit-state function at design point; The Hohenbichler and Rackwitz’s formula also aims at improving the reliability estimate for moderate \(\beta\). It should be noted that the Hohenbichler and Rackwitz’s formula is asymptotically equivalent to Breitung’s formula for a large \(\beta\).

References:

• Hohenbichler, M. and Rackwitz, R., 1988. Improvement of second-order reliability estimates by importance sampling. Journal of Engineering Mechanics, 114(12), pp.2195-2199.

• Hu, Z., Mansour, R., Olsson, M. and Du, X., 2021. Second-order reliability methods: a review and comparative study. Structural and Multidisciplinary Optimization, 64(6), pp.3233-3263.

• Bourinet, J.M., 2018. Reliability analysis and optimal design under uncertainty-Focus on adaptive surrogate-based approaches (Doctoral dissertation, Université Clermont Auvergne).

Function help

from ffpack.rrm import hrackSORM
help( hrackSORM )

Help on function hrackSORM in module ffpack.rrm.secondOrderReliabilityMethod:

hrackSORM(dim, g, dg, distObjs, corrMat, quadDeg=99, quadRange=8, dx=1e-06)
    Second order reliability method based on Hohenbichler and Rackwitz algorithm.

    Parameters
    ----------
    dim: integer
        Space dimension ( number of random variables ).
    g: function
        Limit state function. It will be called like g( [ x1, x2, ... ] ).
    dg: array_like of function
        Gradient of the limit state function. It should be an array_like of function
        like dg = [ dg_dx1, dg_dx2, ... ]. To get the derivative of i-th random
        variable at ( x1*, x2*, ... ), dg[ i ]( x1*, x2*, ... ) will be called.
        dg can be None, see the following Notes.
    distObjs: array_like of distributions
        Marginal distribution objects. It should be the freezed distribution
        objects with pdf, cdf, ppf. We recommend to use scipy.stats functions.
    corrMat: 2d matrix
        Correlation matrix of the marginal distributions.
    quadDeg: integer
        Quadrature degree for Nataf transformation
    quadRange: scalar
        Quadrature range for Nataf transformation. The integral will be performed
        in the range [ -quadRange, quadRange ].
    dx : scalar, optional
        Spacing for auto differentiation. Not required if dg is provided.

    Returns
    -------
    beta: scalar
        Reliability index.
    pf: scalar
        Probability of failure.
    uCoord: 1d array
        Design point coordinate in U space.
    xCoord: 1d array
        Design point coordinate in X space.

    Raises
    ------
    ValueError
        If the dim is less than 1.
        If the dim does not match the disObjs and corrMat.
        If corrMat is not 2d matrix.
        If corrMat is not positive definite.
        If corrMat is not symmetric.
        If corrMat diagonal is not 1.

    Notes
    -----
    If dg is None, the numerical differentiation will be used. The tolerance of the
    numerical differentiation can be changed in globalConfig.

    Examples
    --------
    >>> from ffpack.rrm import tvedtSORM
    >>> dim = 2
    >>> g = lambda X: -np.sum( X ) + 1
    >>> dg = [ lambda X: -1, lambda X: -1 ]
    >>> distObjs = [ stats.norm(), stats.norm() ]
    >>> corrMat = np.eye( dim )
    >>> beta, pf, uCoord, xCoord = hrackSORM( dim, g, dg, distObjs, corrMat )

Example with explicit derivative of LSF

# Define the dimension for the FORM problem
dim = 2

# Define the limit state function (LSF) g
g = lambda X: X[ 0 ] ** 4 + 2 * X[ 1 ] ** 4 - 20

# Explicit derivative of LSF
# dg is a list in which each element is a partial derivative function of g w.r.t. X
# dg[0] = partial g / partial X[0]
# dg[1] = partial g / partial X[1]
dg = [ lambda X: 4 * X[ 0 ] ** 3, lambda X: 8 * X[ 1 ] ** 3 ]

# Marginal distributions and correlation Matrix of the random variables
distObjs = [ stats.norm( 5.0, 5.0 ), stats.norm( 5.0, 5.0 ) ]
corrMat = np.eye( dim )

beta, pf, uCoord, xCoord = hrackSORM( dim, g, dg, distObjs, corrMat )

print( "Reliability index: " )
print( beta )
print()
print( "Failure probability: " )
print( pf )
print()
print( "Design point coordinate in U space: " )
print( uCoord )
print()
print( "Design point coordinate in X space: " )
print( xCoord )

Reliability index:
0.9519628114174661

Failure probability:
0.09883455895117937

Design point coordinate in U space:
[-0.6398897170780924, -0.7048222075811494]

Design point coordinate in X space:
[1.8005514146095387, 1.4758889620942535]

Example with automatic differentiation of LSF

# Define the dimension for the FORM problem
dim = 2

# Define the limit state function (LSF) g
g = lambda X: X[ 0 ] ** 4 + 2 * X[ 1 ] ** 4 - 20

# If dg is None, the internal automatic differentiation algorithm will be used
dg = None

# Marginal distributions and correlation Matrix of the random variables
distObjs = [ stats.norm( 5.0, 5.0 ), stats.norm( 5.0, 5.0 ) ]
corrMat = np.eye( dim )

beta, pf, uCoord, xCoord = hrackSORM( dim, g, dg, distObjs, corrMat )

print( "Reliability index: " )
print( beta )
print()
print( "Failure probability: " )
print( pf )
print()
print( "Design point coordinate in U space: " )
print( uCoord )
print()
print( "Design point coordinate in X space: " )
print( xCoord )

Reliability index:
0.9519628114174661

Failure probability:
0.09883455894748124

Design point coordinate in U space:
[-0.6398897170780924, -0.7048222075811494]

Design point coordinate in X space:
[1.8005514146095387, 1.4758889620942535]