G2 Example 17

G2Example17.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any kind.
# ------------------------------------------------------------------------

from G2AnalysisLinearStatic import *
from G2Model import *


#    ---- F ---> |----------------------- L,E,I --------------------------|
#          ----> |                                                        |
#          ----> |                                                        |
#          ----> |                                                        |
#          ----> |                                                        |
#        q ----> |                                                  H,E,I |
#          ----> |                                                        |
#          ----> |                                                        |
#          ----> |                                                        |
#          ----> |                                                        |
#          ----> |                                                        |
#             -----                                                    -----
#
# -----------------------------------------------------\-------------------
# RANDOM VARIABLES [N, m]
# ------------------------------------------------------------------------
h = 0.25 # Cross-section height and width
#                  E              I              q             F
means = np.array([60e9,          h**4/12,       5000,         15000])
stdvs = np.array([0.15*means[0], 0.05*means[1], 0.2*means[2], 0.3*means[3]])
correlation = []
# ------------------------------------------------------------------------
# LIMIT-STATE FUNCTION
# ------------------------------------------------------------------------
def g(x):

    # Count the total number of evaluations
    global evaluationCount
    evaluationCount += 1

    # Get random variable realizations
    E = x[0]
    I = x[1]
    q = x[2]
    F = x[3]

    # Set value of deterministic variables
    H = 4
    L = 7
    A = h**2

    # Set displacement threshold
    threshold = H/360

    # Nodal coordinates
    NODES = [[0.0, 0.0],
             [0.0, H],
             [L,   H],
             [L,   0.0]]

    # Boundary conditions (0=free, 1=fixed, sets #DOFs per node)
    CONSTRAINTS = [[1, 1, 1],
                   [0, 0, 0],
                   [0, 0, 0],
                   [1, 1, 1]]

    # Element information
    ELEMENTS = [[5, E, A, I, q,   1, 2],
                [5, E, A, I, 0.0, 2, 3],
                [5, E, A, I, 0.0, 3, 4]]

    # Empty arrays
    SECTIONS = np.zeros((3, 3))
    MATERIALS = np.zeros((3, 3))

    # Nodal loads
    LOADS = [[0.0, 0.0, 0.0],
             [F,   0.0, 0.0],
             [0.0, 0.0, 0.0],
             [0.0, 0.0, 0.0]]

    # Mass
    MASS = [[0, 0, 0],
            [0, 0, 0],
            [0, 0, 0],
            [0, 0, 0]]

    # Create the model object
    a = [NODES, CONSTRAINTS, ELEMENTS, SECTIONS, MATERIALS, LOADS, MASS]
    m = model(a)

    # Request response sensitivities calculated with the direct differentiation method (DDM)
    DDMparameters = [['Element', 'E', [1, 2, 3]],
                     ['Element', 'I', [1, 2, 3]],
                     ['Element', 'q', [1]],
                     ['Nodal load', 2, 1]]

    # Analyze
    trackNode = 2
    trackDOF = 1
    u, dudx = linearStaticAnalysis(m, trackNode, trackDOF, DDMparameters)

    # Set the value of the limit-state function and its gradient vector
    g = threshold-u
    dgdx = -np.array(dudx)

    # Issue warning if the limit-state function is negative at the mean: the resulting beta would fool you
    print('\n'"The value of the limit-state function is %.5f" % g)
    if g < 0 and evaluationCount == 1:
        print("   ... and that value is NEGATIVE, implying a negative reliability index, NOT reflected in the results")
        print("       (remove this warning message if you know what you are doing)")
        import sys
        sys.exit()

    return g, dgdx

# ------------------------------------------------------------------------
# ALGORITHM TO FIND THE DESIGN POINT
# ------------------------------------------------------------------------

ddm = True
maxSteps = 50
convergenceCriterion1 = 0.001
convergenceCriterion2 = 0.001
merit_y = 1.0
merit_g = 5.0

# Correlation matrix
numRV = len(means)
R = np.identity(numRV)
for i in range(len(correlation)):
    rv_i = int(correlation[i][0])
    rv_j = int(correlation[i][1])
    rho = correlation[i][2]
    R[rv_i-1, rv_j-1] = rho
    R[rv_j-1, rv_i-1] = rho

# Cholesky decomposition
L = np.linalg.cholesky(R)

# Set start point in the standard normal space
numRV = len(means)
y = np.zeros(numRV)

# Initialize the vectors
x = np.zeros(numRV)
dGdy = np.zeros(numRV)

# Start the FORM loop
plotStepValues = []
plotBetaValues = []
evaluationCount = 0
for i in range(1, maxSteps+1):

    # Transform from y to x
    x = means + np.dot(np.dot(np.diag(stdvs), L), y)

    # Set the Jacobian
    dxdy = np.dot(np.diag(stdvs), L)

    # Evaluate the limit-state function, g(x) = G(y)
    gValue, dgdx = g(x)

    # Gradient in the y-space
    dGdy = dgdx.dot(dxdy)

    # Determine the alpha-vector
    alpha = np.multiply(dGdy, -1 / np.linalg.norm(dGdy))

    # Calculate the first convergence criterion
    if i == 1:
        gfirst = gValue
    criterion1 = np.abs(gValue / gfirst)

    # Calculate the second convergence criterion
    yNormOr1 = np.linalg.norm(y)
    plotStepValues.append(i)
    plotBetaValues.append(yNormOr1)
    if yNormOr1 < 1.0:
        yNormOr1 = 1.0
    u_scaled = np.multiply(y, 1.0/yNormOr1)
    u_scaled_times_alpha = u_scaled.dot(alpha)
    criterion2 = np.linalg.norm(np.subtract(u_scaled, np.multiply(alpha, u_scaled_times_alpha)))

    # Print status
    print('\n'"FORM Step:", i, "Check1:", criterion1, ", Check2:", criterion2)

    # Check convergence
    if criterion1 < convergenceCriterion1 and criterion2 < convergenceCriterion2:

        # Here we converged; first calculate beta and pf
        betaFORM = np.linalg.norm(y)
        print('\n'"FORM analysis converged after %i steps with reliability index %.2f" % (i, betaFORM))
        print('\n'"Total number of limit-state evaluations:", evaluationCount)
        print('\n'"Design point in x-space:", x)
        print('\n'"Design point in y-space:", y)

        # Importance vectors alpha and gamma
        print('\n'"Importance vector alpha:", alpha)
        gamma = alpha.dot(np.linalg.inv(L))
        gamma = gamma / np.linalg.norm(gamma)
        print('\n'"Importance vector gamma:", gamma)

        # Plot the evolution of the analysis
        plt.clf()
        plt.title("Convergence of the Reliability Index")
        plt.grid(True)
        plt.autoscale(True)
        plt.ylabel('Distance from origin to trial point')
        plt.xlabel('Iteration number')
        plt.plot(plotStepValues, plotBetaValues, 'ks-', linewidth=1.0)
        print('\n'"Click somewhere in the plot to continue...")
        plt.waitforbuttonpress()
        break

    # Take a step if convergence did not occur
    else:

        # Give message if non-convergence in maxSteps
        if i == maxSteps:
            print('\n'"FORM analysis did not converge")
            break

        # Determine the search direction
        searchDirection = np.multiply(alpha, (gValue / np.linalg.norm(dGdy) + alpha.dot(y))) - y

        # Define the merit function for the line search along the search direction
        def meritFunction(stepSize):

            # Determine the trial point corresponding to the step size
            yTrial = y + stepSize * searchDirection

            # Calculate the distance from the origin
            yNormTrial = np.linalg.norm(yTrial)

            # Transform from y to x
            xTrial = means + np.multiply(L.dot(y), stdvs)

            # Evaluate the limit-state function
            gTrial, dudx = g(xTrial)

            # Evaluate the merit function m = 1/2 * ||y||^2  + c*g
            return merit_y * yNormTrial + merit_g * np.abs(gTrial/gfirst)

        # Perform the line search for the optimal step size
        stepSize = 1.0

        # Take a step
        y += stepSize * searchDirection

    # Increment the loop counter
    i += 1