.. 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-54 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 than 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 56-58 1. Define some useful functions ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 60-77 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: - with mean :math:`[bt, b(t+\delta t)]` and - with 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 79-92 .. code-block:: Python 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 93-95 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 98-103 .. code-block:: Python 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 104-105 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 108-119 .. code-block:: Python 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 120-132 .. code-block:: Python 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 133-134 This function evaluates the probability using the Low Discrepancy sampling. .. GENERATED FROM PYTHON SOURCE LINES 137-152 .. code-block:: Python 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 153-154 This function evaluates the probability using the FORM algorithm for event systems.. .. GENERATED FROM PYTHON SOURCE LINES 157-164 .. code-block:: Python 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 165-167 2. Evaluate the probability --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 171-176 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 178-199 .. code-block:: Python 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 200-206 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 208-227 .. code-block:: Python 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 228-232 .. code-block:: Python 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.00030517578125, 0.00030517578125, 0.00030517578125, 0.00030517578125, 0.00030517578125, 0.0006103515625, 0.00030517578125, 0.0006103515625, 0.00030517578125, 0.00030517578125, 0.0, 0.0, 0.00030517578125, 0.00030517578125, 0.0, 0.00091552734375, 0.0006103515625, 0.00030517578125, 0.00030517578125, 0.00091552734375, 0.00091552734375, 0.0, 0.00091552734375, 0.0006103515625, 0.00091552734375] 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.407247221976507e-05, 7.202731340749856e-05, 8.087457586090277e-05, 9.070184981054664e-05, 0.00010160352566970394, 0.00011368175073412876, 0.0001270463113912958, 0.00014181491000870818, 0.00015811435561607578, 0.00017607979215241996, 0.00019585595886454746, 0.0002175971126689752, 0.0002414674407523497, 0.0002676410518535999, 0.0002963031343688264, 0.0003276489835870634, 0.00036188514225255717, 0.00039922842206094566, 0.0004399070455775222, 0.000484160922259269, 0.0005322401306923784, 0.0005844062178174964, 0.0006409303353549688, 0.000702094564935671, 0.0007681919118323028, 0.0008395236033398269] .. GENERATED FROM PYTHON SOURCE LINES 233-234 Draw the graphs! .. GENERATED FROM PYTHON SOURCE LINES 236-252 .. code-block:: Python 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("upper left") 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 2.148 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-jupyter :download:`Download Jupyter notebook: plot_proba_system_event.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_proba_system_event.py `