This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # ------------------------------------------------------------------------ | |
| # 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 |