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 r(t) = R - bt and a time-varying load 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:

g(t) \leq 0

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

We propose the following probabilistic model:

  • R is the initial resistance, and R \sim \mathcal{N}(\mu_R, \sigma_R);

  • b is the deterioration rate of the resistance; it is deterministic;

  • S(\omega,t) is the time-varying stress, which is modeled by a stationary Gaussian process of mean value \mu_S, standard deviation \sigma_S and a squared exponential covariance model C(s,t).

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

\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 t, we note the random variable Z_t = R - bt - S_t where S_t = S(., t).

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

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

  • \mathcal{E}_t^1 = \{   Z_t \geq 0\} and

  • \mathcal{E}_t^2 = \{   Z_{t+\Delta t} \leq 0\}.

The objective is to evaluate:

\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 Y_t = (bt + S_t, b(t+\Delta t) + S_{t+\Delta t}). Here, Y_t is a bivariate Normal random vector:

  • whith mean [bt, b(t+\delta t)] and

  • whith covariance matrix \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 Y_t =(Y_t^1, Y_t^2).

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)

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

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])

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

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
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

This function evaluates the probability using the Low Discrepancy sampling.

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

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

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

2. Evaluate the probability

First, fix some parameters: (\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 \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 t and t+\Delta t. We advice to take \Delta t \simeq 10^{-1}.

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

Use all the methods previously described:

  • Monte Carlo: values in values_MC

  • Low discrepancy suites: values in values_QMC

  • FORM: values in values_FORM

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)
    )
print("Values MC = ", values_MC)
print("Values QMC = ", values_QMC)
print("Values FORM = ", values_FORM)
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]

Draw the graphs!

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()
plot proba system event

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