Ticket #146: cantilever_beam_FORMSQPReliability.py

#! /usr/bin/env python

from openturns import *

from math import *

from openturns_viewer import ViewImage


# This function enables a pretty print of a NumericalPoint
def printNumericalPoint(point, digits) :
  oss = "["
  eps = pow(0.1, digits)
  for i in range(point.getDimension()) :
    if i == 0 :
      sep = ""
    else :
      sep = ","
    if fabs(point[i]) < eps :
      oss += sep + str(fabs(point[i]))
    else :
      oss += sep + str(point[i])
    sep = ","
  oss += "]" 
  return oss


###########################################
### Fonction 'poutre'
###########################################

# We create a numerical math function 
stochasticDimension = 4

# Analytical construction : Input
inputFunction = Description(stochasticDimension)
inputFunction[0] = "E"
inputFunction[1] = "F"
inputFunction[2] = "L"
inputFunction[3] = "I"
        
# Analytical construction : Output      
outputFunction = Description(1)
outputFunction[0] = "y"

formulas = Description(outputFunction.getSize())
formulas[0] = "F*(L^3)/(3*E*I)"
        
myFunction = NumericalMathFunction(inputFunction, outputFunction, formulas)



###########################################
### Random input vector
###########################################


dim = myFunction.getInputNumericalPointDimension()

# We create a normal distribution point of dimension 4
mean = NumericalPoint(dim, 0.0)
# E : steel : 210000 MPa
mean[0] = 2.1e11
# F : 1kg : 10N
mean[1] =  10.0
# L : 1 m
mean[2] = 1.0
# I : square hollow section of width 2 cm and thickness 1mm : 2.47325 e-9
mean[3] =  2.47325e-9
sigma = NumericalPoint(dim, 1.0)
# E : 5% * mean
sigma[0] = 0.05 *  mean[0]
# F : 10% * mean
sigma[1] = 0.1 * mean[1]
# L : 1% * mean
sigma[2] = 0.01 * mean[2]
# I : 1% * mean
sigma[3] = 0.01 * mean[3]

R = IdentityMatrix(dim)
myDistribution = Normal(mean, sigma, R)


input = RandomVector(myDistribution)

output =  RandomVector(myFunction, input)

# Custom gradient
eps = NumericalPoint(4, 1.0)
for i in range(4):
        eps[i] = mean[i]*0.001
myGradient = NonCenteredFiniteDifferenceGradient(eps, myFunction.getEvaluationImplementation())
myFunction.setGradientImplementation(NonCenteredFiniteDifferenceGradient(myGradient))


#################################################################
### Probabilistic Study : threshold exceedance: deviation <-1cm
#################################################################

# We create an Event from this RandomVector
threshold = 0.01
myEvent = Event(output, ComparisonOperator(GreaterOrEqual()), threshold)

######
# FORM
######

print "#####"
print "FORM"
print "#####"


# We create a NearestPoint algorithm 
mySQP = SQP()
specific = SQPSpecificParameters(0.5,1.0e-03,1.2)
mySQP.setSpecificParameters(SQPSpecificParameters())
mySQP.setMaximumIterationsNumber(1000)
mySQP.setMaximumAbsoluteError(1.0e-3) # |(ui+1)-ui|
mySQP.setMaximumRelativeError(1.0e-3) # |norm(u)|
mySQP.setMaximumResidualError(1.0e-3) # crit FERUM
mySQP.setMaximumConstraintError(1.0e-3) 
 

# We create a FORM algorithm 
myAlgoFORM = FORM(NearestPointAlgorithm(mySQP), myEvent, mean)


# Perform the simulation 
myAlgoFORM.run()


# Stream out the result 
resultFORM = myAlgoFORM.getResult()
digits = 5
print  "FORM event probability=" , resultFORM.getEventProbability() 
print  "generalized reliability index=" , resultFORM.getGeneralisedReliabilityIndex()