Fault trees and reliability block diagrams

/**
 * Fault trees and reliability block diagrams on one coherent-logic core.
 *   fault tree node:  { event: 'id' }  |  { and: [nodes] }  |  { or: [nodes] }  |  { koon: [nodes], k }   (k or more of the children occur)
 *   events: { id: { p } | { lambda, t } | { lambda, tau } | { mttf, mttr } }
 * The same event id used twice is the SAME event (a repeated event); probability is then computed exactly by factoring (Shannon expansion),
 * so repeated events never get multiplied as if independent. Cut sets: MOCUS (Fussell & Vesely 1972) with absorption.
 */
export function eventProbability(e) {
  if (e == null) return NaN;
  if (e.p != null) return e.p;
  if (e.lambda != null && e.t != null) return 1 - Math.exp(-e.lambda * e.t);                       // non-repairable, failed by time t
  if (e.lambda != null && e.tau != null) { const x = e.lambda * e.tau; return x < 1e-9 ? x / 2 : 1 - (1 - Math.exp(-x)) / x; } // standby tested every tau: mean unavailability
  if (e.mttf != null && e.mttr != null) return e.mttr / (e.mttf + e.mttr);                         // repairable, steady state
  return NaN;
}

/** P(k or more of independent events with probabilities ps occur): the Poisson-binomial tail by dynamic programming. */
export function atLeast(k, ps) {
  let d = [1]; for (const p of ps) { const nd = new Array(d.length + 1).fill(0); for (let i = 0; i < d.length; i++) { nd[i] += d[i] * (1 - p); nd[i + 1] += d[i] * p; } d = nd; }
  let s = 0; for (let i = Math.max(0, k); i < d.length; i++) s += d[i]; return s;
}

function prob(n, p) {
  switch (kind(n)) {
    case 'event': return p.get(n.event);
    case 'and': return kids(n).reduce((a, c) => a * prob(c, p), 1);
    case 'or': return 1 - kids(n).reduce((a, c) => a * (1 - prob(c, p)), 1);
    default: return atLeast(n.k, kids(n).map((c) => prob(c, p)));
  }
}

function exact(n, p) {
  if (n === true) return 1; if (n === false) return 0;
  const ls = leaves(n), count = new Map(); for (const l of ls) count.set(l, (count.get(l) || 0) + 1);
  const rep = [...count].filter(([, c]) => c > 1).map(([id]) => id);
  if (!rep.length) return prob(n, p);
  const x = rep[0], px = p.get(x);
  return px * exact(condition(n, x, true), p) + (1 - px) * exact(condition(n, x, false), p);
}

const csP = cs.sets.map((s) => s.reduce((a, e) => a * p.get(e), 1)), rare = Math.min(1, csP.reduce((a, b) => a + b, 0)), mcub = 1 - csP.reduce((a, b) => a * (1 - b), 1);
if (cs.truncated) probs.push(problem('cutsets-truncated', 'warning', 'There are more cut sets than the limit; the list and its bounds are incomplete.', []));

if (top != null) for (const id of count.keys()) {
  const hi = exact(tree, new Map(p).set(id, 1)), lo = exact(tree, new Map(p).set(id, 0)), pi = p.get(id);
  imp[id] = { p: pi, birnbaum: hi - lo, criticality: P > 0 ? (hi - lo) * pi / P : 0, fussellVesely: P > 0 ? (P - lo) / P : 0, raw: P > 0 ? hi / P : null, rrw: lo > 0 ? P / lo : null };
}

export function compReliability(c, t) {
  if (c.r != null) return c.r;
  if (c.lambda != null) return Math.exp(-c.lambda * t);
  if (c.weibull) return Math.exp(-Math.pow(t / c.weibull.eta, c.weibull.beta));
  return NaN;
}

export function rbd(node, comps, t) {
  const events = {}, probs = [];
  for (const [id, c] of Object.entries(comps)) { const r = compReliability(c, t); events[id] = { p: 1 - r }; if (!(r >= 0 && r <= 1)) probs.push(problem('rbd-r', 'error', `Block "${id}" gives a reliability of ${r}.`, [{ id }])); }
  const koonBad = (n) => (n.koon && !(Number.isInteger(n.k) && n.k >= 1 && n.k <= n.koon.length)) || [n.series, n.parallel, n.koon].some((c) => c && c.some(koonBad));
  if (koonBad(node)) probs.push(problem('rbd-koon', 'error', 'A k-of-n block needs k between 1 and n.', []));
  if (probs.length) return { ok: false, problems: probs, method: M_RBD };
  const ft = faultTree(toFailure(node), events);
  if (!ft.ok) return { ok: false, problems: ft.problems, method: M_RBD };
  return { ok: true, reliability: 1 - ft.top, unreliability: ft.top, cutSets: ft.cutSets, importance: ft.importance, problems: ft.problems, method: M_RBD };
}

/** Mean time to failure of a non-repairable diagram: ∫ R(t) dt by Simpson's rule to the time R falls below 1e-9 (or `horizon`). */
export function rbdMttf(node, comps, { horizon = null, steps = 4000 } = {}) {
  let H = horizon;
  if (H == null) { H = 1; while (H < 1e12 && rbd(node, comps, H).reliability > 1e-9) H *= 2; }
  const n = steps - (steps % 2), h = H / n; let s = 0;
  for (let i = 0; i <= n; i++) { const r = rbd(node, comps, i * h).reliability; s += r * (i === 0 || i === n ? 1 : i % 2 ? 4 : 2); }
  return { mttf: s * h / 3, horizon: H, method: M_RBD };
}

/** Beta-factor common cause (IEC 61508-6 Annex D, simple form): the part of a unit's failure rate shared by the whole redundant group. */
export function betaFactor(lambda, beta) { return { independent: lambda * (1 - beta), common: lambda * beta }; }