Comparing FORM, SORM and simulation#

Define a stochastic model, run FORM, and use its result to compare SORM and simulation estimates.

Before you start: basic Python and NumPy. This example needs only the core PySTRA installation.

Problem and assumptions#

The limit state is

\[g(\gamma,X_1,X_2,X_3)=\gamma-\frac{X_2}{1000X_3}-\left(\frac{X_1}{200X_3}\right)^2.\]

Failure is \(g\leq0\), and \(\gamma=1\). The model inputs are:

Variable	Distribution	Mean	Standard deviation
\(X_1\)	Lognormal	500	100
\(X_2\)	Normal	2000	400
\(X_3\)	Uniform	5	0.5

The specified physical Pearson correlations are 0.3 between \(X_1\) and \(X_2\), and 0.2 for the other pairs. The default Nataf model uses a Gaussian copula.

Define the model#

[1]:

import numpy as np
import pandas as pd

import pystra as ra


def limit_state(gamma, X1, X2, X3):
    return gamma - X2 / (1000 * X3) - (X1 / (200 * X3)) ** 2


model = ra.StochasticModel()
model.add_variable(ra.Lognormal("X1", 500, 100))
model.add_variable(ra.Normal("X2", 2000, 400))
model.add_variable(ra.Uniform("X3", 5, 0.5))
model.add_variable(ra.Constant("gamma", 1))
model.set_correlation(
    ra.CorrelationMatrix(
        [
            [1.0, 0.3, 0.2],
            [0.3, 1.0, 0.2],
            [0.2, 0.2, 1.0],
        ]
    )
)
response = ra.LimitState(limit_state)

Variable names match the limit-state arguments. Distribution constructors use mean and standard deviation by default; they do not use variance. Constants are supplied explicitly through the same model.

Run FORM#

Run the analysis explicitly and inspect its result snapshot. Check convergence before using its probability or design point.

[2]:

form = ra.FORM(model=model, limit_state=response)
result = form.run()
assert result.converged
print(f"Reliability index: {result.beta:.4f}")
print(f"Failure probability: {result.failure_probability:.6g}")
print("Physical design point:", result.design_point_x)
print("Standard-normal design point:", result.design_point_u)

Reliability index: 1.7540
Failure probability: 0.0397173
Physical design point: [ 631.50413468 2310.35249517    4.51737361]
Standard-normal design point: [ 1.27804451  0.40781853 -1.12991998]

FORM approximates the transformed boundary by a tangent plane. Its convergence establishes a local numerical solution; it does not establish the accuracy of that probability approximation.

Compare SORM approximations#

SORM uses local curvature. Supply the converged FORM analysis to reuse the design point. Curve fitting estimates a Hessian; point fitting samples the boundary on each side of the design point.

[3]:

sorm = ra.SORM(model=model, limit_state=response, form=form)
sorm_result = sorm.run()
print(sorm_result.summary())

sorm_points = ra.SORM(model=model, limit_state=response, form=form, options=ra.SORMOptions(fit="point"))
sorm_points_result = sorm_points.run()
print(sorm_points_result.summary())

SORM result
  Status                   converged
  Message                  Fitted
  Failure probability      0.0322193
  Reliability index        1.84913
  Fit                      curve
  FORM reliability index   1.75398
  Pf, breitung             0.0322193
  Pf, modified breitung    0.031148
  Limit-state evaluations  16
SORM result
  Status                   converged
  Message                  Fitted
  Failure probability      0.033365
  Reliability index        1.83349
  Fit                      point
  FORM reliability index   1.75398
  Pf, breitung             0.033365
  Pf, modified breitung    0.0323876
  Limit-state evaluations  52

The reports distinguish Breitung’s approximation and its Hohenbichler–Rackwitz modification. Differences between formulas are approximation differences, not sampling uncertainty.

Check with simulation#

Use a reproducible simulation check of the original event. These existing simulation algorithms use NumPy’s global random state. Their sampling uncertainty must be considered when comparing results.

[4]:

options = ra.SimulationOptions(n_samples=20_000)
monte_carlo = ra.CrudeMonteCarlo(
    model=model,
    limit_state=response,
    options=options,
    rng=2026,
)
monte_carlo_result = monte_carlo.run()
importance = ra.ImportanceSampling(
    model=model,
    limit_state=response,
    options=options,
    rng=2027,
)
importance_result = importance.run()

pd.DataFrame(
    [
        {"Method": "FORM", "Failure probability": result.failure_probability},
        {
            "Method": "Crude Monte Carlo",
            "Failure probability": monte_carlo_result.failure_probability,
        },
        {
            "Method": "Importance sampling",
            "Failure probability": importance_result.failure_probability,
        },
    ]
)

[4]:

	Method	Failure probability
0	FORM	0.039717
1	Crude Monte Carlo	0.032667
2	Importance sampling	0.032037

Interpretation#

Compare the FORM/SORM approximations with the simulation estimates and their reported precision. A close result on this example does not establish accuracy for a different failure geometry.

Continue with SciPy distributions or simulation methods. The FORM and SORM theory explains the geometric assumptions.

Continue: User guide · API reference · Theory