RP110 analysis and 2D graphics

The objective of this example is to present problem 110 of the BBRC. We also present graphic elements for the visualization of the limit state surface in 2 dimensions.

import openturns as ot
import openturns.viewer as otv
import otbenchmark as otb

problem = otb.ReliabilityProblem110()
print(problem)

name = RP110
event = class=ThresholdEventImplementation antecedent=class=CompositeRandomVector function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x0,x1,gsys] evaluationImplementation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0,x1] outputVariablesNames=[gsys] formulas=[if (x0 <= 3.5)
var g1 := 0.85 - 0.1 * x0;
else
  g1 := 4 - x0;
if (x1 <= 2.0)
  var g2 := 2.3 - x1;
else
  g2 := 0.5 - 0.1 * x1;
gsys := min(g1, g2);] gradientImplementation=class=SymbolicGradient name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0,x1] outputVariablesNames=[gsys] formulas=[if (x0 <= 3.5)
var g1 := 0.85 - 0.1 * x0;
else
  g1 := 4 - x0;
if (x1 <= 2.0)
  var g2 := 2.3 - x1;
else
  g2 := 0.5 - 0.1 * x1;
gsys := min(g1, g2);] hessianImplementation=class=SymbolicHessian name=Unnamed evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x0,x1] outputVariablesNames=[gsys] formulas=[if (x0 <= 3.5)
var g1 := 0.85 - 0.1 * x0;
else
  g1 := 4 - x0;
if (x1 <= 2.0)
  var g2 := 2.3 - x1;
else
  g2 := 0.5 - 0.1 * x1;
gsys := min(g1, g2);] antecedent=class=UsualRandomVector distribution=class=JointDistribution name=JointDistribution dimension=2 copula=class=IndependentCopula name=IndependentCopula dimension=2 marginal[0]=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] marginal[1]=class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] operator=class=Less name=Unnamed threshold=0
probability = 3.19e-05

event = problem.getEvent()
g = event.getFunction()

problem.getProbability()

3.19e-05

Create the Monte-Carlo algorithm

algoProb = ot.ProbabilitySimulationAlgorithm(event)
algoProb.setMaximumOuterSampling(1000)
algoProb.setMaximumCoefficientOfVariation(0.01)
algoProb.run()

Get the results

resultAlgo = algoProb.getResult()
neval = g.getEvaluationCallsNumber()
print("Number of function calls = %d" % (neval))
pf = resultAlgo.getProbabilityEstimate()
print("Failure Probability = %.4f" % (pf))
level = 0.95
c95 = resultAlgo.getConfidenceLength(level)
pmin = pf - 0.5 * c95
pmax = pf + 0.5 * c95
print("%.1f %% confidence interval :[%.4f,%.4f] " % (level * 100, pmin, pmax))

Number of function calls = 1000
Failure Probability = 0.0000
95.0 % confidence interval :[0.0000,0.0000]

Compute the bounds of the domain

inputVector = event.getAntecedent()
distribution = inputVector.getDistribution()
X1 = distribution.getMarginal(0)
X2 = distribution.getMarginal(1)
alphaMin = 0.0000001
alphaMax = 1 - alphaMin
lowerBound = ot.Point(
    [X1.computeQuantile(alphaMin)[0], X2.computeQuantile(alphaMin)[0]]
)
upperBound = ot.Point(
    [X1.computeQuantile(alphaMax)[0], X2.computeQuantile(alphaMax)[0]]
)

nbPoints = [100, 100]
figure = g.draw(lowerBound, upperBound, nbPoints)
figure.setTitle(" Iso-values of limit state function")
_ = otv.View(figure)

Print the iso-values of the distribution

_ = otv.View(distribution.drawPDF())

sampleSize = 10000
drawEvent = otb.DrawEvent(event)

cloud = drawEvent.drawSampleCrossCut(sampleSize)
_ = otv.View(cloud)

Draw the limit state surface

bounds = ot.Interval(lowerBound, upperBound)
bounds

class=Interval name=Unnamed dimension=2 lower bound=class=Point name=Unnamed dimension=2 values=[-5.19934,-5.19934] upper bound=class=Point name=Unnamed dimension=2 values=[5.19934,5.19934] finite lower bound=[1,1] finite upper bound=[1,1]

graph = drawEvent.drawLimitStateCrossCut(bounds)
graph.add(cloud)
_ = otv.View(graph)

domain = drawEvent.fillEventCrossCut(bounds)
_ = otv.View(domain)

domain.add(cloud)
_ = otv.View(domain)

otv.View.ShowAll()