Create a new reliability problem

The objective of this example is to create a new reliability problem. This can be useful when we want to consider a new problem which is not already available in the benchmark. This makes it possible to consider a new problem as if it was in the package.

import openturns as ot
import openturns.viewer as otv
import otbenchmark as otb
import numpy as np

In order to create our own benchmark problem, we must create a new class which derives from ReliabilityBenchmarkProblem.

class OscillatorProblem(otb.ReliabilityBenchmarkProblem):
    def __init__(self):
        """Create the Oscillator problem

        This is a two-degree-of-freedom damped oscillator with primary
        and secondary systems.

        Reference
        ---------
        - "Analyse de sensibilité fiabiliste avec prise en compte d'incertitudes
          sur le modèle probabiliste", Vincent Chabridon. Thèse. 2018, pp.91-92.
        - De Stefano M., Der Kiureghian A. (1990). An efficient algorithm for
          second-order reliability analysis.
          Report No. UCB/SEMM-90/20. Dept of Civil and Environmental Engineering,
          University of California, Berkeley.
        """

        def oscillator(x):
            mp, ms, kp, ks, xip, xis, S0 = x
            omega_p = np.sqrt(kp / mp)
            omega_s = np.sqrt(ks / ms)
            gamma = ms / mp
            omega_a = 0.5 * (omega_p + omega_s)
            xi_a = 0.5 * (xip + xis)
            theta = 1.0 / omega_a * (omega_p - omega_s)
            factor_1 = np.pi * S0 / (4.0 * xis * omega_s**3)
            factor_2 = (
                xi_a * xis / (xip * xis * (4.0 * xi_a**2 + theta**2) + gamma * xi_a**2)
            )
            factor_3 = (
                (xip * omega_p**3 + xis * omega_s**3)
                * omega_p
                / (4.0 * xi_a * omega_a**4)
            )
            F = 3 * ks * np.sqrt(factor_1 * factor_2 * factor_3)
            return [F]

        name = "Oscillator"
        dim = 7
        limitStateFunction = ot.PythonFunction(dim, 1, oscillator)
        mean_list = [1.5, 0.01, 1.0, 0.01, 0.05, 0.02, 100.0]
        cov_list = [0.1, 0.1, 0.2, 0.2, 0.4, 0.5, 0.1]
        myCollection = ot.DistributionCollection(dim)
        for i, (mu, cov) in enumerate(zip(mean_list, cov_list)):
            parameters = ot.LogNormalMuSigma(mu, mu * cov, 0.0)
            myCollection[i] = ot.ParametrizedDistribution(parameters)
        distribution = ot.ComposedDistribution(myCollection)
        inputRandomVector = ot.RandomVector(distribution)
        outputRandomVector = ot.CompositeRandomVector(
            limitStateFunction, inputRandomVector
        )
        threshold = 0.0
        thresholdEvent = ot.ThresholdEvent(outputRandomVector, ot.Less(), threshold)
        probability = 3.78e-7
        super().__init__(name, thresholdEvent, probability)
        return None

problem = OscillatorProblem()
print(problem)

name = Oscillator
event = class=ThresholdEventImplementation antecedent=class=CompositeRandomVector function=class=Function name=Unnamed implementation=class=FunctionImplementation name=Unnamed description=[x0,x1,x2,x3,x4,x5,x6,y0] evaluationImplementation=class=PythonEvaluation name=OpenTURNSPythonFunction description=[x0,x1,x2,x3,x4,x5,x6,y0] parameter=class=Point name=Unnamed dimension=0 values=[] gradientImplementation=class=CenteredFiniteDifferenceGradient name=Unnamed epsilon=class=Point name=Unnamed dimension=7 values=[1e-05,1e-05,1e-05,1e-05,1e-05,1e-05,1e-05] evaluation=class=PythonEvaluation name=OpenTURNSPythonFunction description=[x0,x1,x2,x3,x4,x5,x6,y0] parameter=class=Point name=Unnamed dimension=0 values=[] hessianImplementation=class=CenteredFiniteDifferenceHessian name=Unnamed epsilon=class=Point name=Unnamed dimension=7 values=[0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001] evaluation=class=PythonEvaluation name=OpenTURNSPythonFunction description=[x0,x1,x2,x3,x4,x5,x6,y0] parameter=class=Point name=Unnamed dimension=0 values=[] antecedent=class=UsualRandomVector distribution=class=JointDistribution name=JointDistribution dimension=7 copula=class=IndependentCopula name=IndependentCopula dimension=7 marginal[0]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=1.5 sigma=0.15 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=0.40049 sigmaLog=0.0997513 gamma=0 marginal[1]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=0.01 sigma=0.001 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=-4.61015 sigmaLog=0.0997513 gamma=0 marginal[2]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=1 sigma=0.2 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=-0.0196104 sigmaLog=0.198042 gamma=0 marginal[3]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=0.01 sigma=0.002 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=-4.62478 sigmaLog=0.198042 gamma=0 marginal[4]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=0.05 sigma=0.02 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=-3.06994 sigmaLog=0.385253 gamma=0 marginal[5]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=0.02 sigma=0.01 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=-4.02359 sigmaLog=0.472381 gamma=0 marginal[6]=class=ParametrizedDistribution parameters=class=LogNormalMuSigma name=Unnamed mu=100 sigma=10 gamma=0 distribution=class=LogNormal name=LogNormal dimension=1 muLog=4.6002 sigmaLog=0.0997513 gamma=0 operator=class=Less name=Unnamed threshold=0
probability = 3.78e-07

event = problem.getEvent()
g = event.getFunction()

problem.getProbability()

3.78e-07

Create the Monte-Carlo algorithm

algoProb = ot.ProbabilitySimulationAlgorithm(event)
algoProb.setMaximumOuterSampling(100)
algoProb.setBlockSize(10)
algoProb.setMaximumCoefficientOfVariation(0.01)
algoProb.run()

Get the results

resultAlgo = algoProb.getResult()
neval = g.getEvaluationCallsNumber()
print("Number of function calls = %d" % (neval))
pf = resultAlgo.getProbabilityEstimate()
print("Failure Probability = %.4f" % (pf))
level = 0.95
c95 = resultAlgo.getConfidenceLength(level)
pmin = pf - 0.5 * c95
pmax = pf + 0.5 * c95
print("%.1f %% confidence interval :[%.4f,%.4f] " % (level * 100, pmin, pmax))

Number of function calls = 1000
Failure Probability = 0.0000
95.0 % confidence interval :[0.0000,0.0000]

The reference probability is too close to zero: Monte-Carlo sampling cannot achieve a satisfactory accuracy.

otv.View.ShowAll()