.. 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 Click :ref:here  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 from openturns import * 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) .. 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 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! .. GENERATED FROM PYTHON SOURCE LINES 95-103 .. 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) .. 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 107-123 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 124-125 This function evaluates the probability using the Low Discrepancy sampling. .. GENERATED FROM PYTHON SOURCE LINES 127-144 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 145-146 This function evaluates the probability using the FORM algorithm for event systems.. .. GENERATED FROM PYTHON SOURCE LINES 148-161 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 162-164 2. Evaluate the probability --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 168-172 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}. .. GENERATED FROM PYTHON SOURCE LINES 174-195 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 196-202 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 204-216 .. 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 217-221 .. 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.00030517578125, 0.0, 0.0, 0.0, 0.0, 0.0, 0.00030517578125, 0.00030517578125, 0.0, 0.00091552734375, 0.0006103515625, 0.00030517578125, 0.00091552734375, 0.001220703125, 0.00030517578125, 0.001220703125, 0.00091552734375, 0.0006103515625, 0.00030517578125, 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.407247221452685e-05, 7.202731340860951e-05, 8.087457491593016e-05, 9.070179169300293e-05, 0.0001016035263802752, 0.00011368175169084065, 0.00012704623305305574, 0.00014181490835112135, 0.00015811426182631293, 0.00017607968850372457, 0.00019585584543730799, 0.00021759698560570485, 0.0002414672698574692, 0.00026764105252706364, 0.0002963031350828803, 0.0003276489830651007, 0.00036188490016252284, 0.00039922815388919713, 0.00043990704675780126, 0.00048416092659680056, 0.0005322401297909951, 0.0005844058510196042, 0.0006409299329991489, 0.0007020945699336272, 0.0007681919182910387, 0.0008395236089949951] .. GENERATED FROM PYTHON SOURCE LINES 222-223 Draw the graphs! .. GENERATED FROM PYTHON SOURCE LINES 225-241 .. 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") 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 8.768 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_proba_system_event.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_proba_system_event.py  .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: plot_proba_system_event.ipynb  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery _