Estimate a flooding probability

In this example, we estimate the probability that the output of a function exceeds a given threshold with the FORM method. We consider the flooding model.

Define the model

from openturns.usecases import flood_model
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt

ot.Log.Show(ot.Log.NONE)

We load the flooding model from the usecases module :

fm = flood_model.FloodModel()
distribution = fm.distribution
model = fm.model.getMarginal(1)

See the input distribution

distribution

JointDistribution

• name=JointDistribution
• dimension: 8
• description=[Q (m3/s),Ks,Zv (m),Zm (m),B (m),L (m),Zb (m),Hd (m)]
• copula: IndependentCopula(dimension = 8)

Index	Variable	Distribution
0	Q (m3/s)	TruncatedDistribution(Gumbel(beta = 558, gamma = 1013), bounds = [0, (19000.8) +inf[)
1	Ks	TruncatedDistribution(Normal(mu = 30, sigma = 7.5), bounds = [0, (87.3797) +inf[)
2	Zv (m)	Uniform(a = 49, b = 51)
3	Zm (m)	Uniform(a = 54, b = 56)
4	B (m)	Triangular(a = 295, m = 300, b = 305)
5	L (m)	Triangular(a = 4990, m = 5000, b = 5010)
6	Zb (m)	Triangular(a = 55, m = 55.5, b = 56)
7	Hd (m)	Uniform(a = 2, b = 4)

See the model

model.getOutputDescription()

[S]

Draw the distribution of a sample of the output.

sampleSize = 1000
inputSample = distribution.getSample(sampleSize)
outputSample = model(inputSample)
graph = ot.HistogramFactory().build(outputSample).drawPDF()
_ = viewer.View(graph)

Define the event

Then we create the event whose probability we want to estimate.

vect = ot.RandomVector(distribution)
G = ot.CompositeRandomVector(model, vect)
event = ot.ThresholdEvent(G, ot.Greater(), 0.0)
event.setName("overflow")

Estimate the probability with FORM

Define a solver.

optimAlgo = ot.Cobyla()
optimAlgo.setMaximumCallsNumber(1000)
optimAlgo.setMaximumAbsoluteError(1.0e-8)
optimAlgo.setMaximumRelativeError(1.0e-10)
optimAlgo.setMaximumResidualError(1.0e-10)
optimAlgo.setMaximumConstraintError(1.0e-10)

Run FORM.

startingPoint = distribution.getMean()
algo = ot.FORM(optimAlgo, event, startingPoint)
algo.run()
result = algo.getResult()
standardSpaceDesignPoint = result.getStandardSpaceDesignPoint()

Retrieve results.

result = algo.getResult()
probability = result.getEventProbability()
print("Pf=", probability)

Pf= 0.000650134456773035

Importance factors.

graph = result.drawImportanceFactors()
view = viewer.View(graph)
plt.show()