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

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}\]

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)

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

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

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

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)

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

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)

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.

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

(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}\).

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)

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