.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_auto_reliability_sensitivity_reliability_plot_probaSystemEvent.py>`     to download the full example code
    .. rst-class:: sphx-glr-example-title

    .. _sphx_glr_auto_reliability_sensitivity_reliability_plot_probaSystemEvent.py:


Time variant system reliability problem
=======================================

The objective is to evaluate the outcrossing rate from a safe to a failure domain in a time variant reliability problem.

We consider the following limit state function, defined as the difference between a degrading resistance :math:`r(t) = R - bt`  and a time-varying load :math:`S(t)`:

..math:
  \begin{align*}
  g(t)= r(t) - S(t) = R - bt - S(t) \quad \forall t \in [0,T]
  \end{align*}

The failure domaine is defined by:

.. math::
   g(t) \leq 0


which means that the resistance at :math:`t` is less thant the stress at :math:`t`.


We propose the following probabilistic model: 

- :math:`R` is the initial resistance, and :math:`R \sim \mathcal{N}(\mu_R, \sigma_R)`;   
- :math:`b` is the deterioration rate of the resistance; it is deterministic; 
- :math:`S(\omega,t)` is the time-varying stress, which is modeled by a stationary Gaussian process of mean value :math:`\mu_S`, standard deviation :math:`\sigma_S` and a squared exponential covariance model :math:`C(s,t)`.


The outcrossing rate from the safe to the failure domain at instant :math:`t` is defined by:

.. math::
   \nu^+(t) = \lim_{\Delta t \rightarrow 0+} \dfrac{\mathbb{P}\{ g(t) \ge 0 \cap g(t+\Delta t) \leq 0\} }{\Delta t} 


For each :math:`t`, we note the random variable  :math:`Z_t = R - bt - S_t` where :math:`S_t = S(., t)`. 

To evaluate :math:`\nu^+(t)`, we need to consider the bivariate random vector :math:`(Z_t, Z_{t+\Delta t})`.

The event :math:`\{ g(t) \geq 0 \cap g(t+\Delta t) <0\}` writes as the intersection of both events :

- :math:`\mathcal{E}_t^1 = \{   Z_t \geq 0\}` and
- :math:`\mathcal{E}_t^2 = \{   Z_{t+\Delta t} \leq 0\}`.

The objective is to evaluate:

.. math::
   \mathbb{P}\{\mathcal{E}_t^1 \cap \mathcal{E}_t^2\} \quad \forall t \in [0,T]


1. Define some useful functions
-------------------------------

We define the bivariate random vector :math:`Y_t = (bt + S_t, b(t+\Delta t) + S_{t+\Delta t})`. 
Here, :math:`Y_t` is a bivariate Normal random vector: 

- whith mean :math:`[bt, b(t+\delta t)]` and 
- whith covariance matrix :math:`\Sigma` defined by:

..math::
  \begin{align*}
  \Sigma = \left(
  \begin{array}{cc}
  C(t, t) & C(t, t+\Delta t) \\
  C(t, t+\Delta t) & C(t+\Delta t, t+\Delta t)
  \end{array}
  \right)
  \end{align*}

This function buils :math:`Y_t =(Y_t^1, Y_t^2)`.


.. code-block:: default

    def buildNormal(b, t, mu_S, covariance, delta_t = 1e-5):
        sigma = CovarianceMatrix(2)
        sigma[0, 0] = covariance(t, t)[0,0]
        sigma[0, 1] = covariance(t, t+delta_t)[0,0]
        sigma[1, 1] = covariance(t+delta_t, t+delta_t)[0,0]
        return Normal([b*t + mu_S, b*(t+delta_t) + mu_S], sigma)









This function creates the trivariate random vector :math:`(R, Y_t^1, Y_t^2)` where :math:`R` is independant from :math:`(Y_t^1, Y_t^2)`. We need to create this random vector because both events  :math:`\mathcal{E}_t^1` and :math:`\mathcal{E}_t^2` must be defined from the same random vector!


.. code-block:: default

    def buildCrossing(b, t, mu_S, covariance, R, delta_t = 1e-5):
        normal = buildNormal(b, t, mu_S, covariance, delta_t)
        #return BlockIndependentDistribution([R, normal]): only from the 1.16 version!
        marg = [R, normal.getMarginal(0), normal.getMarginal(1)]
        cop = ComposedCopula([IndependentCopula(1), normal.getCopula()])
        return ComposedDistribution(marg, cop)









This function evaluates the probability using the Monte Carlo sampling. It defines the intersection event :math:`\mathcal{E}_t^1 \cap \mathcal{E}_t^2`.


.. code-block:: default

    def computeCrossingProbability_MonteCarlo(b, t, mu_S, covariance, R, delta_t, n_block, n_iter, CoV):
        full = buildCrossing(b, t, mu_S, covariance, R, delta_t)
        X = RandomVector(full)
        f1 = SymbolicFunction(["R", "X1", "X2"], ["X1 - R"])
        e1 = ThresholdEvent(CompositeRandomVector(f1, X), Less(), 0.0)
        f2 = SymbolicFunction(["R", "X1", "X2"], ["X2 - R"])
        e2 = ThresholdEvent(CompositeRandomVector(f2, X), GreaterOrEqual(), 0.0)
        event = IntersectionEvent([e1, e2])
        algo = ProbabilitySimulationAlgorithm(event, MonteCarloExperiment())
        algo.setBlockSize(n_block)
        algo.setMaximumOuterSampling(n_iter)
        algo.setMaximumCoefficientOfVariation(CoV)
        algo.run()
        return algo.getResult().getProbabilityEstimate() / delta_t









This function evaluates the probability using the Low Discrepancy sampling. 


.. code-block:: default

    def computeCrossingProbability_QMC(b, t, mu_S, covariance, R, delta_t, n_block, n_iter, CoV):
        full = buildCrossing(b, t, mu_S, covariance, R, delta_t)
        X = RandomVector(full)
        f1 = SymbolicFunction(["R", "X1", "X2"], ["X1 - R"])
        e1 = ThresholdEvent(CompositeRandomVector(f1, X), Less(), 0.0)
        f2 = SymbolicFunction(["R", "X1", "X2"], ["X2 - R"])
        e2 = ThresholdEvent(CompositeRandomVector(f2, X), GreaterOrEqual(), 0.0)
        event = IntersectionEvent([e1, e2])
        algo = ProbabilitySimulationAlgorithm(event, LowDiscrepancyExperiment(SobolSequence(X.getDimension()), n_block, False))
        algo.setBlockSize(n_block)
        algo.setMaximumOuterSampling(n_iter)
        algo.setMaximumCoefficientOfVariation(CoV)
        algo.run()
        return algo.getResult().getProbabilityEstimate() / delta_t









This function evaluates the probability using the FORM algorithm for event systems..


.. code-block:: default

    def computeCrossingProbability_FORM(b, t, mu_S, covariance, R, delta_t):
        full = buildCrossing(b, t, mu_S, covariance, R, delta_t)
        X = RandomVector(full)
        f1 = SymbolicFunction(["R", "X1", "X2"], ["X1 - R"])
        e1 = ThresholdEvent(CompositeRandomVector(f1, X), Less(), 0.0)
        f2 = SymbolicFunction(["R", "X1", "X2"], ["X2 - R"])
        e2 = ThresholdEvent(CompositeRandomVector(f2, X), GreaterOrEqual(), 0.0)
        event = IntersectionEvent([e1, e2])
        algo = SystemFORM(SQP(), event, X.getMean())
        algo.run()
        return algo.getResult().getEventProbability() / delta_t









2. Evaluate the probability
---------------------------


.. code-block:: default

    from openturns import *
    from openturns.viewer import View
    from math import sqrt








First, fix some parameters: :math:`(\mu_R, \sigma_R, \mu_S, \sigma_S, \Delta t, T, b)` and the covariance model wich is the Squared Exponential model.
Be careful to the parameter  :math:`\Delta t` which is of great importance: if it is too small, the simulation methods have problems to converge because the correlation rate is too high between the instants :math:`t` and :math:`t+\Delta t`.
We advice to take :math:`\Delta t \simeq 10^{-1}`.



.. code-block:: default

    mu_S = 3.0
    sigma_S = 0.5
    l = 10

    b = 0.01

    mu_R = 5.0
    sigma_R = 0.3
    R = Normal(mu_R, sigma_R)

    covariance = SquaredExponential([l/sqrt(2)], [sigma_S])

    t0 = 0.0
    t1 = 50.0
    N = 26

    # Get all the time steps t
    times = RegularGrid(t0, (t1 - t0) / (N - 1.0), N).getVertices()

    delta_t = 1e-1








Use all the methods previously described:

- Monte Carlo: values in values_MC
- Low discrepancy suites: values in values_QMC
- FORM: values in values_FORM



.. code-block:: default

    values_MC = list()
    values_QMC = list()
    values_FORM = list()

    for tick in times:
        values_MC.append(computeCrossingProbability_MonteCarlo(b, tick[0], mu_S, covariance, R, delta_t, 2**12, 2**3, 1e-2))
        values_QMC.append(computeCrossingProbability_QMC(b, tick[0], mu_S, covariance, R, delta_t, 2**12, 2**3, 1e-2))
        values_FORM.append(computeCrossingProbability_FORM(b, tick[0], mu_S, covariance, R, delta_t))









.. code-block:: default

    print('Values MC = ', values_MC)
    print('Values QMC = ', values_QMC)
    print('Values FORM = ', values_FORM)





.. rst-class:: sphx-glr-script-out

 Out:

 .. code-block:: none

    Values MC =  [0.0, 0.0, 0.0, 0.00030517578125, 0.0, 0.0, 0.00030517578125, 0.0, 0.00030517578125, 0.0, 0.00030517578125, 0.0, 0.00030517578125, 0.00030517578125, 0.0006103515625, 0.0006103515625, 0.0, 0.0, 0.0006103515625, 0.00091552734375, 0.00030517578125, 0.0, 0.001220703125, 0.00030517578125, 0.0006103515625, 0.0006103515625]
    Values QMC =  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0006103515625, 0.0, 0.0006103515625, 0.0, 0.0, 0.0, 0.00030517578125, 0.00030517578125, 0.0, 0.00030517578125, 0.0, 0.00091552734375, 0.0006103515625, 0.0, 0.00030517578125, 0.0, 0.0006103515625, 0.001220703125, 0.00091552734375, 0.0006103515625]
    Values FORM =  [6.407247221452685e-05, 7.202731340860951e-05, 8.087457491593016e-05, 9.070179169300293e-05, 0.0001016035263802752, 0.00011368175169091608, 0.00012704623305297141, 0.00014181490835112135, 0.00015811426182631293, 0.00017607968850349097, 0.0001958558454373012, 0.00021759698560569734, 0.0002414672698574692, 0.00026764105252706364, 0.0002963031350828803, 0.0003276489830651007, 0.00036188490016252284, 0.00039922815388919713, 0.0004399070467586194, 0.00048416092659680056, 0.0005322401297909951, 0.0005844058510196042, 0.0006409299329987239, 0.0007020945699345352, 0.0007681919182910387, 0.0008395236089949951]




Draw the graphs!


.. code-block:: default

    g = Graph()
    g.setAxes(True)
    g.setGrid(True)
    c = Curve(times, [[p] for p in values_MC])
    g.add(c)
    c = Curve(times, [[p] for p in values_QMC])
    g.add(c)
    c = Curve(times, [[p] for p in values_FORM])
    g.add(c)
    g.setLegends(["MC", "QMC", "FORM"])
    g.setColors(["red", "blue", 'black'])
    g.setLegendPosition("topleft")
    g.setXTitle("t")
    g.setYTitle("Outcrossing rate")
    Show(g)



.. image:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_probaSystemEvent_001.png
    :alt: plot probaSystemEvent
    :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  5.927 seconds)


.. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_probaSystemEvent.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_probaSystemEvent.py <plot_probaSystemEvent.py>`



  .. container:: sphx-glr-download sphx-glr-download-jupyter

     :download:`Download Jupyter notebook: plot_probaSystemEvent.ipynb <plot_probaSystemEvent.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_