.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here ` to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_reliability_sensitivity_reliability_plot_event_system.py:
System events: unions or intersections of events
================================================
This example illustrates system events, which are defined as unions or intersections of other events. We will show how to estimate their probability both with Monte-Carlo sampling (using the class `ProbabilitySimulationAlgorithm`) and with a first order approximation (using the class `SystemFORM`).
**Intersection**
The event defined as the intersection of several events is realized when all sub-events occurs:
.. math::
E_{sys} = \bigcap_{i=1}^N E_i
**Union**
The event defined as the union of several events is realized when at least one sub-event occurs:
.. math::
E_{sys} = \bigcup_{i=1}^N E_i
.. code-block:: default
from __future__ import print_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
For system events, we always have to use the same root cause.
.. code-block:: default
dim = 2
distribution = ot.Normal(dim)
X = ot.RandomVector(distribution)
Define some basic events E1, E2 and E3.
.. code-block:: default
f1 = ot.SymbolicFunction(['x0', 'x1'], ['x0'])
f2 = ot.SymbolicFunction(['x0', 'x1'], ['x1'])
f3 = ot.SymbolicFunction(['x0', 'x1'], ['x0+x1'])
Y1 = ot.CompositeRandomVector(f1, X)
Y2 = ot.CompositeRandomVector(f2, X)
Y3 = ot.CompositeRandomVector(f3, X)
e1 = ot.ThresholdEvent(Y1, ot.Less(), 0.0) # E1 <=> x0<0
e2 = ot.ThresholdEvent(Y2, ot.Greater(), 0.0) # E2 <=> x1>0
e3 = ot.ThresholdEvent(Y3, ot.Greater(), 0.0) # E3 <=> x0+x1>0
Define the intersection E3=E1 AND E2.
.. code-block:: default
e4 = ot.IntersectionEvent([e1, e2])
Approximate probability of that event: :math:`\approx 1/4`.
.. code-block:: default
e4.getSample(10000).computeMean()
.. raw:: html
[0.253]
Define the union E4=E1 OR E2.
.. code-block:: default
e5 = ot.UnionEvent([e1, e2])
Approximate probability of that event: :math:`\approx ~3/4`.
.. code-block:: default
e5.getSample(10000).computeMean()
.. raw:: html
[0.7472]
It supports recursion: define E7=E1 OR (E1 AND E3).
.. code-block:: default
e7 = ot.UnionEvent([e1, ot.IntersectionEvent([e2, e3])])
print(e3.isComposite())
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
True
With Monte-Carlo algorithm
--------------------------
Of course, we can use simulation algorithms.
.. code-block:: default
experiment = ot.MonteCarloExperiment()
algo = ot.ProbabilitySimulationAlgorithm(e7, experiment)
algo.setMaximumOuterSampling(2500)
algo.setBlockSize(4)
algo.setMaximumCoefficientOfVariation(-1.0)
algo.run()
result = algo.getResult()
result.getProbabilityEstimate()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
0.7479999999999998
SystemFORM
----------
`SystemFORM` is an approximation method suitable for system events.
The event must be in its disjunctive normal form (union of intersections, or a single intersection).
root cause
.. code-block:: default
dim = 5
mean = [200.0] * dim
mean[-1] = 60
mean[-2] = 60
sigma = [30.0] * dim
sigma[-1] = 15.0
R = ot.CorrelationMatrix(dim)
for i in range(dim):
for j in range(i):
R[i, j] = 0.5
dist = ot.Normal(mean, sigma, R)
leaf events
.. code-block:: default
X = ot.RandomVector(dist)
inputs = ['M1', 'M2', 'M3', 'M4', 'M5']
e0 = ot.ThresholdEvent(ot.CompositeRandomVector(ot.SymbolicFunction(inputs, ['M1-M2+M4']), X), ot.Less(), 0.0)
e1 = ot.ThresholdEvent(ot.CompositeRandomVector(ot.SymbolicFunction(inputs, ['M2+2*M3-M4']), X), ot.Less(), 0.0)
e2 = ot.ThresholdEvent(ot.CompositeRandomVector(ot.SymbolicFunction(inputs, ['2*M3-2*M4-M5']), X), ot.Less(), 0.0)
e3 = ot.ThresholdEvent(ot.CompositeRandomVector(ot.SymbolicFunction(inputs, ['-(M1+M2+M4+M5-5*10.0)']), X), ot.Less(), 0.0)
e4 = ot.ThresholdEvent(ot.CompositeRandomVector(ot.SymbolicFunction(inputs, ['-(M2+2*M3+M4-5*40.0)']), X), ot.Less(), 0.0)
system event in DNF form (union of intersections)
.. code-block:: default
event = ot.UnionEvent([ot.IntersectionEvent([e0, e3, e4]), ot.IntersectionEvent([e2, e3, e4])])
compute probability with basic sampling
.. code-block:: default
event.getSample(10000).computeMean()
.. raw:: html
[0.0828]
Run systemFORM.
.. code-block:: default
solver = ot.AbdoRackwitz()
solver.setMaximumIterationNumber(1000)
solver.setMaximumAbsoluteError(1.0e-3)
solver.setMaximumRelativeError(1.0e-3)
solver.setMaximumResidualError(1.0e-3)
solver.setMaximumConstraintError(1.0e-3)
algo = ot.SystemFORM(solver, event, mean)
algo.run()
result = algo.getResult()
result.getEventProbability()
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
0.07883551213326333
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 0.246 seconds)
.. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_event_system.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_event_system.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_event_system.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_