.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_proba_system_event.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_reliability_sensitivity_reliability_plot_proba_system_event.py: Time variant system reliability problem ======================================= .. GENERATED FROM PYTHON SOURCE LINES 8-53 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] .. GENERATED FROM PYTHON SOURCE LINES 55-57 1. Define some useful functions ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 59-76 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)`. .. GENERATED FROM PYTHON SOURCE LINES 78-91 .. code-block:: default from math import sqrt from openturns.viewer import View import openturns as ot def buildNormal(b, t, mu_S, covariance, delta_t=1e-5): sigma = ot.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 ot.Normal([b * t + mu_S, b * (t + delta_t) + mu_S], sigma) .. GENERATED FROM PYTHON SOURCE LINES 92-93 This function creates the trivariate random vector :math:`(R, Y_t^1, Y_t^2)` where :math:`R` is independent 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! .. GENERATED FROM PYTHON SOURCE LINES 95-100 .. 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 ot.BlockIndependentDistribution([R, normal]) .. GENERATED FROM PYTHON SOURCE LINES 101-102 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`. .. GENERATED FROM PYTHON SOURCE LINES 104-115 .. code-block:: default def getXEvent(b, t, mu_S, covariance, R, delta_t): full = buildCrossing(b, t, mu_S, covariance, R, delta_t) X = ot.RandomVector(full) f1 = ot.SymbolicFunction(["R", "X1", "X2"], ["X1 - R"]) e1 = ot.ThresholdEvent(ot.CompositeRandomVector(f1, X), ot.Less(), 0.0) f2 = ot.SymbolicFunction(["R", "X1", "X2"], ["X2 - R"]) e2 = ot.ThresholdEvent(ot.CompositeRandomVector(f2, X), ot.GreaterOrEqual(), 0.0) event = ot.IntersectionEvent([e1, e2]) return X, event .. GENERATED FROM PYTHON SOURCE LINES 116-128 .. code-block:: default def computeCrossingProbability_MonteCarlo( b, t, mu_S, covariance, R, delta_t, n_block, n_iter, CoV ): X, event = getXEvent(b, t, mu_S, covariance, R, delta_t) algo = ot.ProbabilitySimulationAlgorithm(event, ot.MonteCarloExperiment()) algo.setBlockSize(n_block) algo.setMaximumOuterSampling(n_iter) algo.setMaximumCoefficientOfVariation(CoV) algo.run() return algo.getResult().getProbabilityEstimate() / delta_t .. GENERATED FROM PYTHON SOURCE LINES 129-130 This function evaluates the probability using the Low Discrepancy sampling. .. GENERATED FROM PYTHON SOURCE LINES 132-147 .. code-block:: default def computeCrossingProbability_QMC( b, t, mu_S, covariance, R, delta_t, n_block, n_iter, CoV ): X, event = getXEvent(b, t, mu_S, covariance, R, delta_t) algo = ot.ProbabilitySimulationAlgorithm( event, ot.LowDiscrepancyExperiment(ot.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 .. GENERATED FROM PYTHON SOURCE LINES 148-149 This function evaluates the probability using the FORM algorithm for event systems.. .. GENERATED FROM PYTHON SOURCE LINES 151-158 .. code-block:: default def computeCrossingProbability_FORM(b, t, mu_S, covariance, R, delta_t): X, event = getXEvent(b, t, mu_S, covariance, R, delta_t) algo = ot.SystemFORM(ot.SQP(), event, X.getMean()) algo.run() return algo.getResult().getEventProbability() / delta_t .. GENERATED FROM PYTHON SOURCE LINES 159-161 2. Evaluate the probability --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 165-169 First, fix some parameters: :math:`(\mu_R, \sigma_R, \mu_S, \sigma_S, \Delta t, T, b)` and the covariance model which 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}`. .. GENERATED FROM PYTHON SOURCE LINES 171-192 .. code-block:: default mu_S = 3.0 sigma_S = 0.5 ll = 10 b = 0.01 mu_R = 5.0 sigma_R = 0.3 R = ot.Normal(mu_R, sigma_R) covariance = ot.SquaredExponential([ll / sqrt(2)], [sigma_S]) t0 = 0.0 t1 = 50.0 N = 26 # Get all the time steps t times = ot.RegularGrid(t0, (t1 - t0) / (N - 1.0), N).getVertices() delta_t = 1e-1 .. GENERATED FROM PYTHON SOURCE LINES 193-199 Use all the methods previously described: - Monte Carlo: values in values_MC - Low discrepancy suites: values in values_QMC - FORM: values in values_FORM .. GENERATED FROM PYTHON SOURCE LINES 201-220 .. 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) ) .. GENERATED FROM PYTHON SOURCE LINES 221-225 .. code-block:: default print("Values MC = ", values_MC) print("Values QMC = ", values_QMC) print("Values FORM = ", values_FORM) .. rst-class:: sphx-glr-script-out .. code-block:: none Values MC = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00030517578125, 0.0, 0.0006103515625, 0.0, 0.0, 0.0, 0.0006103515625, 0.0, 0.0006103515625, 0.0006103515625, 0.00030517578125, 0.00030517578125, 0.0, 0.0006103515625, 0.0006103515625, 0.00030517578125, 0.001220703125, 0.00091552734375, 0.0006103515625, 0.001220703125] 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.407247215635151e-05, 7.202731335077214e-05, 8.087457564322269e-05, 9.070185003762527e-05, 0.00010160352566792341, 0.00011368175043642132, 0.0001270463113567539, 0.00014181490999748653, 0.00015811435561321604, 0.00017607979141239492, 0.00019585595855828544, 0.0002175971122863016, 0.0002414674411439194, 0.0002676410529682008, 0.0002963031348934401, 0.0003276489827287258, 0.0003618851406896481, 0.0003992284220557085, 0.00043990704747756264, 0.0004841609222492271, 0.0005322401306591526, 0.0005844062178788379, 0.0006409303353549256, 0.0007020945630748503, 0.0007681919142532408, 0.0008395236026959429] .. GENERATED FROM PYTHON SOURCE LINES 226-227 Draw the graphs! .. GENERATED FROM PYTHON SOURCE LINES 229-245 .. code-block:: default g = ot.Graph() g.setAxes(True) g.setGrid(True) c = ot.Curve(times, [[p] for p in values_MC]) g.add(c) c = ot.Curve(times, [[p] for p in values_QMC]) g.add(c) c = ot.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") view = View(g) view.ShowAll() .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_proba_system_event_001.png :alt: plot proba system event :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_proba_system_event_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 1.787 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_proba_system_event.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_proba_system_event.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_proba_system_event.ipynb `