System Reliability#

This notebook demonstrates the pystra.system topology layer on classical system-reliability examples. The key distinction is that Pystra asks the user to encode the system event they intend to analyse; it does not infer the structural topology automatically.

The examples cover:

Example

Reference

Purpose

Four-branch case

Schueremans2005

Series system as a minimum of branch limit states

Daniels bundle

Daniels1945

Classical equal-load-sharing benchmark

k-of-n system

Ditlevsen2007

Event-counting topology

Cut-set system

Song2003

Minimal cut-set encoding

Ditlevsen bounds and Mainçon series benchmark

Ditlevsen1979, Maincon2000

Bounds from component and pairwise event probabilities

 
[1]:

import itertools
import sys
from pathlib import Path

import numpy as np
from scipy.stats import binom

for candidate in (Path.cwd(), *Path.cwd().parents):
    src = candidate / "src"
    if (src / "pystra").exists():
        sys.path.insert(0, str(src))
        break

import pystra as ra

rng = np.random.default_rng(20260425)

Four-Branch Series System#

The four-branch case is a standard benchmark for algorithms that must find several disconnected failure regions. It is naturally a series system:

\[g_\mathrm{sys}(x) = \min(g_1, g_2, g_3, g_4).\]

For independent \(X_1, X_2 \sim \mathcal{N}(0, 1)\), the branch functions are:

\[\begin{split}\begin{aligned} g_1 &= 3 + 0.1(X_1 - X_2)^2 - \frac{X_1 + X_2}{\sqrt{2}}, \\ g_2 &= 3 + 0.1(X_1 - X_2)^2 + \frac{X_1 + X_2}{\sqrt{2}}, \\ g_3 &= (X_1 - X_2) + \frac{6}{\sqrt{2}}, \\ g_4 &= (X_2 - X_1) + \frac{6}{\sqrt{2}}. \end{aligned}\end{split}\]
[2]:

sqrt2 = np.sqrt(2.0)
four_branch = ra.SeriesSystem(
    [
        ra.Component("g1", lambda X1, X2: 3.0 + 0.1 * (X1 - X2) ** 2 - (X1 + X2) / sqrt2),
        ra.Component("g2", lambda X1, X2: 3.0 + 0.1 * (X1 - X2) ** 2 + (X1 + X2) / sqrt2),
        ra.Component("g3", lambda X1, X2: (X1 - X2) + 6.0 / sqrt2),
        ra.Component("g4", lambda X1, X2: (X2 - X1) + 6.0 / sqrt2),
    ],
    name="four_branch",
)

n_samples = 200_000
x1 = rng.standard_normal(n_samples)
x2 = rng.standard_normal(n_samples)
failed = four_branch.failure_mask(X1=x1, X2=x2)

pf_hat = failed.mean()
cov_hat = np.sqrt((1.0 - pf_hat) / (n_samples * pf_hat))
float(pf_hat), float(cov_hat)

[2]:

(0.00434, 0.0338684769766651)

The estimate should be close to the usual reference probability of about \(4.5 \times 10^{-3}\). The precise Monte Carlo value changes with the sample size and seed; the point here is that the system topology is encoded directly from the branch functions.

Daniels Equal-Load-Sharing Bundle#

Daniels’ six-thread example gives the individual thread strengths

\[3.2,\ 5.2,\ 4.4,\ 8.2,\ 6.1,\ 5.7.\]

For an equal-load-sharing bundle, sort the strengths in descending order and compute \(r\sigma_r\). The bundle strength is the maximum of those products.

[3]:

thread_strengths = np.array([3.2, 5.2, 4.4, 8.2, 6.1, 5.7])
ordered_strengths = np.sort(thread_strengths)[::-1]
remaining_threads = np.arange(1, len(ordered_strengths) + 1)
bundle_loads = remaining_threads * ordered_strengths

np.column_stack([remaining_threads, ordered_strengths, bundle_loads]), float(bundle_loads.max())

[3]:

(array([[ 1. ,  8.2,  8.2],
        [ 2. ,  6.1, 12.2],
        [ 3. ,  5.7, 17.1],
        [ 4. ,  5.2, 20.8],
        [ 5. ,  4.4, 22. ],
        [ 6. ,  3.2, 19.2]]),
 22.0)

The maximum is 22.0, matching Daniels’ published example. This is not a KofNSystem: after each thread breaks, the load redistributes, so the limit state depends on order statistics rather than only a fixed number of failed components.

k-of-n Topology#

A k-of-n system fails when at least k components have failed. This event-counting representation is useful for simulation and exact Boolean enumeration, but it is discontinuous and should not be treated as a smooth FORM/SORM limit state.

[4]:

n_components = 6
k_failures = 3
component_pf = 0.05

components = [
    ra.Component(f"E{i}", lambda event=f"E{i}", **kwargs: kwargs[event])
    for i in range(n_components)
]
k_of_n = ra.KofNSystem(components, k=k_failures, name="three_of_six")

states = np.array(list(itertools.product([False, True], repeat=n_components))).T
state_values = {f"E{i}": np.where(states[i], -1.0, 1.0) for i in range(n_components)}
weights = np.prod(np.where(states, component_pf, 1.0 - component_pf), axis=0)

pf_enumerated = weights[k_of_n.failure_mask(**state_values)].sum()
pf_binomial = binom.sf(k_failures - 1, n_components, component_pf)
float(pf_enumerated), float(pf_binomial)

[4]:

(0.00222984375, 0.00222984375)

The exact enumeration agrees with the binomial result because the component events were independent and identically distributed in this example.

Minimal Cut Sets#

Song and Der Kiureghian’s rigid-plastic cantilever-bar example uses the system failure event

\[E_\mathrm{sys} = E_1E_2 \cup E_3E_4 \cup E_3E_5.\]

Pystra can encode this topology directly. The probability calculation in the paper uses component and joint event probabilities with linear-programming bounds; those LP bounds are a later extension. The topology itself is already usable for simulation when the component limit states are available.

[5]:

cut_components = {
    name: ra.Component(name, lambda event=name, **kwargs: kwargs[event])
    for name in ("E1", "E2", "E3", "E4", "E5")
}

cantilever_system = ra.CutSetSystem(
    [["E1", "E2"], ["E3", "E4"], ["E3", "E5"]],
    components=cut_components,
    name="cantilever_cut_sets",
)

example_events = {
    "E1": np.array([False, True, False, False]),
    "E2": np.array([False, True, False, True]),
    "E3": np.array([True, False, True, True]),
    "E4": np.array([False, False, True, False]),
    "E5": np.array([False, False, False, True]),
}
example_values = {name: np.where(mask, -1.0, 1.0) for name, mask in example_events.items()}

cantilever_system.failure_mask(**example_values)

[5]:

array([False,  True,  True,  True])

The four samples above fail only when at least one of the supplied cut sets is complete.

Ditlevsen Bounds#

Ditlevsen bounds use component event probabilities and pairwise intersections to bound a union of failure events. For two events the result collapses to exact inclusion-exclusion.

[6]:

ra.ditlevsen_bounds([0.1, 0.2], {(0, 1): 0.03})

[6]:

(0.27, 0.27)

Mainçon’s identical 100-element series-system example reports component failure probability \(P_1 = 1.969 \times 10^{-3}\) and pairwise intersection probability \(P_{12} = 1.361 \times 10^{-3}\). With equal pairwise intersections, the Ditlevsen upper bound is approximately \(6.216 \times 10^{-2}\).

[7]:

n_events = 100
p1 = 1.969e-3
p12 = 1.361e-3

probabilities = np.full(n_events, p1)
intersections = np.full((n_events, n_events), p12)
np.fill_diagonal(intersections, p1)

ra.ditlevsen_bounds(probabilities, intersections)

[7]:

(0.0025769999999999994, 0.06216099999999983)

These bounds are deliberately separate from the topology classes. A topology tells Pystra which samples fail; Ditlevsen bounds require marginal and pairwise event probabilities, usually obtained from component analyses, simulation, or reference data.

Component-based system FORM and original-system Monte Carlo#

SystemFORM runs each component separately in the same complete stochastic model. It supports series and parallel events, including homogeneous nesting. Mixed and cut/tie-set systems use original-system simulation for now.

The four-branch tangent-plane probability is about 0.00539, whereas the nonlinear event probability is about 0.00446. The difference is FORM approximation error, not sampling error. Bounds reported here apply to the tangent-plane events only.

[8]:

model = ra.StochasticModel()
model.addVariable(ra.Normal("X1", 0.0, 1.0))
model.addVariable(ra.Normal("X2", 0.0, 1.0))

system_form = ra.SystemFORM(four_branch, model)
system_form.run()
print("System FORM:", system_form.getFailure())
print("Linearized-event bounds:", system_form.bounds)
print("Original nonlinear event (earlier simulation):", pf_hat)
print("Component betas:", system_form.betas)

System FORM: 0.005392305084089622
Linearized-event bounds: (0.005392305084089622, 0.005394124771863334)
Original nonlinear event (earlier simulation): 0.00434
Component betas: [3.00042432 2.9995758  3.         3.        ]

For an exact check, two independent components with beta = 2 give series probability 2*p - p*p and parallel probability p*p. Monte Carlo below evaluates the original series event. Its CoV reports sampling uncertainty; the sample budget is a maximum and the run may stop early at the target CoV. No observed failures would require more sampling, not a claim of zero risk.

[9]:

from scipy.stats import norm

linear_model = ra.StochasticModel()
linear_model.addVariable(ra.Normal("X", 0.0, 1.0))
linear_model.addVariable(ra.Normal("Y", 0.0, 1.0))
linear_components = [ra.Component("a", lambda X: 2.0-X),
                     ra.Component("b", lambda Y: 2.0-Y)]
series = ra.SeriesSystem(linear_components)
p = norm.sf(2.0)
for topology, exact in [(series, 2*p-p*p),
                        (ra.ParallelSystem(linear_components), p*p)]:
    result = ra.SystemFORM(topology, linear_model)
    result.run()
    print(topology.name, result.getFailure(), "exact:", exact)
    np.testing.assert_allclose(result.getFailure(), exact, rtol=1e-6)

options = ra.AnalysisOptions()
options.setSamples(20_000)
options.setPrintOutput(False)
np.random.seed(2026)
mc = ra.CrudeMonteCarlo(stochastic_model=linear_model,
                       limit_state=series.as_limit_state(),
                       analysis_options=options)
mc.run()
print("Monte Carlo:", mc.getFailure(), "CoV:", mc.cov_q_bar[mc.k-1],
      "samples:", mc.k)

series 0.04498269539267562 exact: 0.04498269539269883
parallel 0.0005175685036590233 exact: 0.0005175685036595636
Monte Carlo: 0.04288888888888889 CoV: 0.04979520523121099 samples: 9000