Benchmark the reliability solvers on the problems

In this example, we show how to run all the methods on all the problems and get the computed probability.

import openturns as ot
import numpy as np
import otbenchmark as otb
import pandas as pd
from tqdm import tqdm
ot.Log.Show(ot.Log.NONE)

We import the list of reliability problems.

benchmarkProblemList = otb.ReliabilityBenchmarkProblemList()
numberOfProblems = len(benchmarkProblemList)
numberOfProblems
26

For each problem in the benchmark, print the problem name and the exact failure probability.

for i in range(numberOfProblems):
    problem = benchmarkProblemList[i]
    name = problem.getName()
    pf = problem.getProbability()
    print("#", i, " : ", name, ", exact PF : ", pf)
# 0  :  RP8 , exact PF :  0.0007897927545597477
# 1  :  RP14 , exact PF :  0.00077285
# 2  :  RP22 , exact PF :  0.004207305511299618
# 3  :  RP24 , exact PF :  0.00286
# 4  :  RP25 , exact PF :  4.148566293759747e-05
# 5  :  RP28 , exact PF :  1.4532945550025393e-07
# 6  :  RP31 , exact PF :  0.003226681209587691
# 7  :  RP33 , exact PF :  0.00257
# 8  :  RP35 , exact PF :  0.00347894632
# 9  :  RP38 , exact PF :  0.0081
# 10  :  RP53 , exact PF :  0.0313
# 11  :  RP55 , exact PF :  0.5600144282863704
# 12  :  RP54 , exact PF :  0.000998
# 13  :  RP57 , exact PF :  0.0284
# 14  :  RP75 , exact PF :  0.00981929872154689
# 15  :  RP89 , exact PF :  0.00543
# 16  :  RP107 , exact PF :  2.92e-07
# 17  :  RP110 , exact PF :  3.19e-05
# 18  :  RP111 , exact PF :  7.65e-07
# 19  :  RP63 , exact PF :  0.000379
# 20  :  RP91 , exact PF :  0.000697
# 21  :  RP60 , exact PF :  0.0456
# 22  :  RP77 , exact PF :  2.87e-07
# 23  :  Four-branch serial system , exact PF :  0.0022227950661944398
# 24  :  R-S , exact PF :  0.07864960352514257
# 25  :  Axial stressed beam , exact PF :  0.02919819462483095

Run several algorithms on a single problem

We want to run several algorithms on a single problem. We set the parameters of the algorithms and run them on a single problem.

maximumEvaluationNumber = 1000
maximumAbsoluteError = 1.0e-3
maximumRelativeError = 1.0e-3
maximumResidualError = 1.0e-3
maximumConstraintError = 1.0e-3
nearestPointAlgorithm = ot.AbdoRackwitz()
nearestPointAlgorithm.setMaximumCallsNumber(maximumEvaluationNumber)
nearestPointAlgorithm.setMaximumAbsoluteError(maximumAbsoluteError)
nearestPointAlgorithm.setMaximumRelativeError(maximumRelativeError)
nearestPointAlgorithm.setMaximumResidualError(maximumResidualError)
nearestPointAlgorithm.setMaximumConstraintError(maximumConstraintError)
i = 3
problem = benchmarkProblemList[i]
metaAlgorithm = otb.ReliabilityBenchmarkMetaAlgorithm(problem)

We try the FORM algorithm.

benchmarkFORM = metaAlgorithm.runFORM(nearestPointAlgorithm)
s1 = benchmarkFORM.summary()
print(s1)
computedProbability = 0.006209245091320793
exactProbability = 0.00286
absoluteError = 0.003349245091320793
numberOfCorrectDigits = 0.0
numberOfFunctionEvaluations = 6
numberOfDigitsPerEvaluation = 0.0

Then the SORM algorithm.

benchmarkSORM = metaAlgorithm.runSORM(nearestPointAlgorithm)
s2 = benchmarkSORM.summary()
print(s2)
computedProbability = 0.006209245091320793
exactProbability = 0.00286
absoluteError = 0.003349245091320793
numberOfCorrectDigits = 0.0
numberOfFunctionEvaluations = 6
numberOfDigitsPerEvaluation = 0.0
benchmarkMC = metaAlgorithm.runMonteCarlo(
    maximumOuterSampling=1000000, coefficientOfVariation=0.1, blockSize=1,
)
s3 = benchmarkMC.summary()
print(s3)
computedProbability = 0.0032433834976647656
exactProbability = 0.00286
absoluteError = 0.00038338349766476546
numberOfCorrectDigits = 0.8727326179291324
numberOfFunctionEvaluations = 30832
numberOfDigitsPerEvaluation = 2.8306065708651153e-05
benchmarkFORMIS = metaAlgorithm.runFORMImportanceSampling(
    nearestPointAlgorithm,
    maximumOuterSampling=1000,
    coefficientOfVariation=0.1,
    blockSize=1,
)
s4 = benchmarkFORMIS.summary()
print(s4)
computedProbability = 0.0023766831479580456
exactProbability = 0.00286
absoluteError = 0.0004833168520419545
numberOfCorrectDigits = 0.7721340949708474
numberOfFunctionEvaluations = 765
numberOfDigitsPerEvaluation = 0.001009325614340977
benchmarkSS = metaAlgorithm.runSubsetSampling(
    maximumOuterSampling=5000, coefficientOfVariation=0.1, blockSize=1,
)
s5 = benchmarkSS.summary()
print(s5)
computedProbability = 0.0028699999999999997
exactProbability = 0.00286
absoluteError = 9.999999999999593e-06
numberOfCorrectDigits = 2.4563660331290604
numberOfFunctionEvaluations = 15000
numberOfDigitsPerEvaluation = 0.00016375773554193736

Run all algorithms on all problems and produce a single result table

For several algorithms and all the reliability problems, we want to estimate the failure probability and compare them.

We create a list of problem names.

problem_names = []
for i in range(numberOfProblems):
    problem = benchmarkProblemList[i]
    name = problem.getName()
    problem_names.append(name)
metrics = [
    "Exact",
    "FORM",
    "SORM",
    "Monte Carlo",
    "FORM-IS",
    "Subset",
]
results = np.zeros((numberOfProblems, len(metrics)))
maximumOuterSampling = 10 ** 2
blockSize = 10 ** 2
coefficientOfVariation = 0.0

for i in tqdm(range(numberOfProblems)):
    problem = benchmarkProblemList[i]
    results[i][0] = problem.getProbability()
    metaAlgorithm = otb.ReliabilityBenchmarkMetaAlgorithm(problem)
    benchmarkResult = metaAlgorithm.runFORM(nearestPointAlgorithm)
    results[i][1] = benchmarkResult.computedProbability
    benchmarkResult = metaAlgorithm.runSORM(nearestPointAlgorithm)
    results[i][2] = benchmarkResult.computedProbability
    benchmarkResult = metaAlgorithm.runMonteCarlo(
        maximumOuterSampling=maximumOuterSampling,
        coefficientOfVariation=coefficientOfVariation,
        blockSize=blockSize,
    )
    results[i][3] = benchmarkResult.computedProbability
    benchmarkResult = metaAlgorithm.runFORMImportanceSampling(
        nearestPointAlgorithm,
        maximumOuterSampling=maximumOuterSampling,
        coefficientOfVariation=coefficientOfVariation,
        blockSize=blockSize,
    )
    results[i][4] = benchmarkResult.computedProbability
    benchmarkResult = metaAlgorithm.runSubsetSampling(
        maximumOuterSampling=maximumOuterSampling,
        coefficientOfVariation=coefficientOfVariation,
        blockSize=blockSize,
    )
    results[i][5] = benchmarkResult.computedProbability

df = pd.DataFrame(results, index=problem_names, columns=metrics)
# df.to_csv("reliability_benchmark_table-output.csv")
df
  0%|          | 0/26 [00:00<?, ?it/s]
  4%|▍         | 1/26 [00:01<00:42,  1.68s/it]
  8%|▊         | 2/26 [00:02<00:22,  1.07it/s]
 15%|█▌        | 4/26 [00:02<00:08,  2.59it/s]
 23%|██▎       | 6/26 [00:02<00:04,  4.05it/s]
 31%|███       | 8/26 [00:02<00:03,  5.76it/s]
 38%|███▊      | 10/26 [00:02<00:02,  6.70it/s]
 50%|█████     | 13/26 [00:03<00:01,  8.09it/s]
 62%|██████▏   | 16/26 [00:03<00:00, 10.63it/s]
 69%|██████▉   | 18/26 [00:03<00:00,  8.65it/s]
 77%|███████▋  | 20/26 [00:05<00:02,  2.64it/s]
 81%|████████  | 21/26 [00:05<00:01,  2.96it/s]
 85%|████████▍ | 22/26 [00:06<00:01,  2.28it/s]
 88%|████████▊ | 23/26 [00:06<00:01,  2.70it/s]
100%|██████████| 26/26 [00:06<00:00,  4.65it/s]
100%|██████████| 26/26 [00:06<00:00,  3.78it/s]
Exact FORM SORM Monte Carlo FORM-IS Subset
RP8 7.897928e-04 6.598878e-04 7.838036e-04 0.0005 8.194590e-04 8.431000e-04
RP14 7.728500e-04 7.003011e-04 6.995436e-04 0.0007 7.618629e-04 7.775000e-04
RP22 4.207306e-03 6.209672e-03 4.390902e-03 0.0037 4.239900e-03 4.567000e-03
RP24 2.860000e-03 6.209245e-03 6.209245e-03 0.0038 2.791203e-03 2.711000e-03
RP25 4.148566e-05 0.000000e+00 0.000000e+00 0.0000 0.000000e+00 3.877234e-05
RP28 1.453295e-07 2.850470e-08 0.000000e+00 0.0000 1.858184e-07 2.398380e-07
RP31 3.226681e-03 2.275013e-02 2.275013e-02 0.0023 2.939076e-03 3.487000e-03
RP33 2.570000e-03 1.349898e-03 1.349898e-03 0.0023 2.461629e-03 2.451000e-03
RP35 3.478946e-03 1.349898e-03 2.134376e-03 0.0033 2.472850e-03 3.948000e-03
RP38 8.100000e-03 7.902212e-03 8.029356e-03 0.0080 8.246963e-03 7.937000e-03
RP53 3.130000e-02 1.180398e-01 2.986164e-02 0.0313 3.069200e-02 3.337000e-02
RP55 5.600144e-01 0.000000e+00 0.000000e+00 0.5659 0.000000e+00 5.635000e-01
RP54 9.980000e-04 5.555704e-02 3.554811e-03 0.0012 1.013974e-03 9.510480e-04
RP57 2.840000e-02 0.000000e+00 0.000000e+00 0.0283 0.000000e+00 2.748000e-02
RP75 9.819299e-03 0.000000e+00 0.000000e+00 0.0096 0.000000e+00 9.592000e-03
RP89 5.430000e-03 2.008594e-09 2.008594e-09 0.0055 3.011492e-04 5.935000e-03
RP107 2.920000e-07 2.866516e-07 2.866516e-07 0.0000 2.880225e-07 3.467086e-07
RP110 3.190000e-05 3.167124e-05 3.167124e-05 0.0000 3.123988e-05 1.497501e-05
RP111 7.650000e-07 0.000000e+00 0.000000e+00 0.0000 0.000000e+00 8.682309e-07
RP63 3.790000e-04 9.999966e-01 0.000000e+00 0.0003 0.000000e+00 3.889000e-04
RP91 6.970000e-04 6.994296e-04 7.011592e-04 0.0006 7.117930e-04 6.269000e-04
RP60 4.560000e-02 4.483968e-02 4.483968e-02 0.0452 4.575646e-02 4.447000e-02
RP77 2.870000e-07 0.000000e+00 0.000000e+00 0.0000 0.000000e+00 2.706695e-07
Four-branch serial system 2.222795e-03 0.000000e+00 0.000000e+00 0.0020 0.000000e+00 2.020000e-03
R-S 7.864960e-02 7.864960e-02 7.864960e-02 0.0768 7.661601e-02 8.100000e-02
Axial stressed beam 2.919819e-02 2.998280e-02 2.933256e-02 0.0298 2.902224e-02 2.856000e-02


Run several algorithms on all problems and get detailed statistics

Run several algorithms on all reliability benchmark problems: print statistics on each problem.

def FormatRow(benchmarkResult):
    """Format a single row of the benchmark table"""
    result = [
        benchmarkResult.exactProbability,
        benchmarkResult.computedProbability,
        benchmarkResult.absoluteError,
        benchmarkResult.numberOfCorrectDigits,
        benchmarkResult.numberOfFunctionEvaluations,
        benchmarkResult.numberOfDigitsPerEvaluation,
    ]
    return result
method_names = ["Monte-Carlo", "FORM", "SORM", "FORM-IS", "SUBSET"]

maximumOuterSampling = 10 ** 2
blockSize = 10 ** 2
coefficientOfVariation = 0.0

result = dict()
for i in range(numberOfProblems):
    problem = benchmarkProblemList[i]
    name = problem_names[i]
    exact_pf_name = "%10s" % ("Exact PF " + name[0:10])
    metrics = [
        exact_pf_name,
        "Estimated PF",
        "Absolute Error",
        "Correct Digits",
        "Function Calls",
        "Digits / Evaluation",
    ]
    results = np.zeros((len(method_names), len(metrics)))
    metaAlgorithm = otb.ReliabilityBenchmarkMetaAlgorithm(problem)
    # Monte-Carlo
    benchmarkResult = metaAlgorithm.runMonteCarlo(
        maximumOuterSampling=maximumOuterSampling,
        coefficientOfVariation=coefficientOfVariation,
        blockSize=blockSize,
    )
    results[0, :] = FormatRow(benchmarkResult)
    # FORM
    benchmarkResult = metaAlgorithm.runFORM(nearestPointAlgorithm)
    results[1, :] = FormatRow(benchmarkResult)
    # SORM
    benchmarkResult = metaAlgorithm.runSORM(nearestPointAlgorithm)
    results[2, :] = FormatRow(benchmarkResult)
    # FORM-IS
    benchmarkResult = metaAlgorithm.runFORMImportanceSampling(
        nearestPointAlgorithm,
        maximumOuterSampling=maximumOuterSampling,
        coefficientOfVariation=coefficientOfVariation,
        blockSize=blockSize,
    )
    results[3, :] = FormatRow(benchmarkResult)
    # Subset
    benchmarkResult = metaAlgorithm.runSubsetSampling(
        maximumOuterSampling=maximumOuterSampling,
        coefficientOfVariation=coefficientOfVariation,
        blockSize=blockSize,
    )
    results[4, :] = FormatRow(benchmarkResult)
    # Gather statistics and print them
    df = pd.DataFrame(results, index=method_names, columns=metrics,)
    # Format the columns for readability
    s = df.style.format(
        {
            exact_pf_name: lambda x: "{:.3e}".format(x),
            "Estimated PF": lambda x: "{:.3e}".format(x),
            "Absolute Error": lambda x: "{:.3e}".format(x),
            "Correct Digits": lambda x: "{:.1f}".format(x),
            "Function Calls": lambda x: "{:d}".format(int(x)),
            "Digits / Evaluation": lambda x: "{:.1f}".format(x),
        }
    )
    result[name] = s
result["RP33"]
  Exact PF RP33 Estimated PF Absolute Error Correct Digits Function Calls Digits / Evaluation
Monte-Carlo 2.570e-03 1.900e-03 6.700e-04 0.6 10000 0.0
FORM 2.570e-03 1.350e-03 1.220e-03 0.3 26 0.0
SORM 2.570e-03 1.350e-03 1.220e-03 0.3 51 0.0
FORM-IS 2.570e-03 3.441e-03 8.709e-04 0.5 10026 0.0
SUBSET 2.570e-03 2.614e-03 4.400e-05 1.8 30000 0.0


result["RP35"]
  Exact PF RP35 Estimated PF Absolute Error Correct Digits Function Calls Digits / Evaluation
Monte-Carlo 3.479e-03 3.800e-03 3.211e-04 1.0 10000 0.0
FORM 3.479e-03 1.350e-03 2.129e-03 0.2 20 0.0
SORM 3.479e-03 2.134e-03 1.345e-03 0.4 33 0.0
FORM-IS 3.479e-03 2.415e-03 1.064e-03 0.5 10020 0.0
SUBSET 3.479e-03 2.939e-03 5.399e-04 0.8 30000 0.0


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