# SPDX-License-Identifier: MIT """YS-EN-0004 ยท Spares for a one-way trip. How many spares does a fleet of 225 Model 4077 cores need to keep working for 40 years with no resupply? Every number below that describes the 4077 (MTBF, how the failure rate splits across the four modules and the chassis, the relative masses) is an INPUT, not a measurement. Change INPUTS and re-run. Model (constant rates first, then relaxed): item i (modules M1..M4 and the chassis) has n = 225 installed copies, rate lam_i, mission T. With replacement on failure the fleet's demand for item i over the mission is D_i ~ Poisson(n lam_i T) (Poisson additivity; failed parts are replaced by fresh spares). Probability of sufficiency (NASA's POS): POS_i = P(D_i <= s_i). With cannibalisation (a dead core donates its good parts), the number of working cores at the end is W = min(n, min_i (n + s_i - D_i)) so "at least k of n cores working" holds iff D_i <= s_i + (n - k) for every i: tolerating m = n - k dead cores is worth m free spares of every item. Independent items: P(W >= k) = prod_i P(D_i <= s_i + m). Strategies compared, all sized to the same fleet-level target: A spare whole cores, no repair below core level: any failure costs a core, D ~ Poisson(n T / MTBF). B spare whole cores, with cannibalisation: a spare core is a bundle of one of each item, so s cores give s of every item. C spare modules and chassis separately, allocated greedily for least mass (Sherbrooke's marginal analysis). Usage: python sparing.py (writes results.json; requires numpy, scipy) """ from __future__ import annotations import json import math from pathlib import Path import numpy as np from scipy import stats HERE = Path(__file__).resolve().parent INPUTS = { "fleet": 225, # cores aboard (9 crates x 25) "mission_yr": 40.0, "core_mtbf_yr": 40.0, # datasheet figure, whole core # share of the core's failure rate carried by each replaceable item "rate_share": {"M1": 0.30, "M2": 0.25, "M3": 0.20, "M4": 0.20, "chassis": 0.05}, # mass of each item as a fraction of one complete core "mass": {"M1": 0.12, "M2": 0.12, "M3": 0.12, "M4": 0.12, "chassis": 0.52}, "targets": [0.9, 0.99, 0.999], } ITEMS = list(INPUTS["rate_share"]) def demand_means(inp=INPUTS) -> dict[str, float]: """Expected fleet failures of each item over the mission.""" core_rate = 1.0 / inp["core_mtbf_yr"] return {i: inp["fleet"] * core_rate * sh * inp["mission_yr"] for i, sh in inp["rate_share"].items()} def poisson_quantile(mean: float, p: float) -> int: """Smallest s with P(Poisson(mean) <= s) >= p.""" s = int(stats.poisson.ppf(p, mean)) while stats.poisson.cdf(s, mean) < p: # guard against ppf rounding s += 1 return s def greedy_allocation(cdfs: dict[str, callable], mass: dict[str, float], target: float, slack: int = 0) -> dict[str, int]: """Least-mass spares with prod_i cdf_i(s_i + slack) >= target. Marginal analysis: start from zero and repeatedly add the spare that buys the most log-POS per unit mass. Each log-cdf is concave in s for Poisson distributions, so for them greedy is optimal on the efficient frontier; the final step can overshoot by one item. For the mixed (Poisson-lognormal) demand that concavity is not guaranteed and the result is a good allocation rather than a proven optimum. """ s = {i: 0 for i in cdfs} logp = {i: math.log(max(cdfs[i](slack), 1e-300)) for i in cdfs} log_target = math.log(target) while sum(logp.values()) < log_target: best, best_gain, best_new = None, -1.0, None for i in cdfs: new = math.log(max(cdfs[i](s[i] + 1 + slack), 1e-300)) gain = (new - logp[i]) / mass[i] if gain > best_gain: best, best_gain, best_new = i, gain, new s[best] += 1 logp[best] = best_new return s def pos_system(cdfs, s, slack=0) -> float: return math.prod(cdfs[i](s[i] + slack) for i in cdfs) def strategies(inp=INPUTS) -> dict: mu = demand_means(inp) mass = inp["mass"] cdfs = {i: (lambda k, m=mu[i]: stats.poisson.cdf(k, m)) for i in ITEMS} total_mu = sum(mu.values()) out = {"demand_mean": mu, "total_demand_mean": total_mu, "by_target": {}} for p in inp["targets"]: a = poisson_quantile(total_mu, p) # A: every failure costs a core # B: smallest s with prod_i P(D_i <= s) >= p b = 0 while pos_system(cdfs, {i: b for i in ITEMS}) < p: b += 1 c = greedy_allocation(cdfs, mass, p) c_mass = sum(c[i] * mass[i] for i in ITEMS) out["by_target"][str(p)] = { "A_spare_cores": a, "A_mass": float(a), "B_spare_cores": b, "B_mass": float(b), "C_spares": c, "C_mass": c_mass, "C_count": sum(c.values()), "C_pos": pos_system(cdfs, c), } return out def k_of_n(inp=INPUTS, p=0.999, tolerances=(0, 5, 10, 20)) -> dict: """Module spares needed when the fleet may lose m cores and still succeed.""" mu = demand_means(inp) cdfs = {i: (lambda k, m=mu[i]: stats.poisson.cdf(k, m)) for i in ITEMS} rows = {} for m in tolerances: s = greedy_allocation(cdfs, inp["mass"], p, slack=m) rows[str(m)] = {"spares": s, "mass": sum(s[i] * inp["mass"][i] for i in ITEMS), "count": sum(s.values())} return rows # --- rate uncertainty --------------------------------------------------------- def lognormal_sigma(error_factor: float) -> float: """NASA/ISS 'error factor' = 95th / 50th percentile of a lognormal rate.""" return math.log(error_factor) / stats.norm.ppf(0.95) def rate_uncertainty(inp=INPUTS, p=0.999, efs=(1.0, 1.5, 2.0, 3.0, 5.0), draws=20000, seed=4077) -> dict: """Design at the nominal rates, then ask what POS the manifest really has. True rate_i = nominal_i * X_i, X_i lognormal with median 1 and the given error factor. Two cases: 'independent' (each item's error drawn separately) and 'common' (one shared multiplier: the whole design is better or worse than predicted, as with a single test campaign). Also sizes the manifest two hedged ways: 'predictive' POS computed against the mixed (Poisson-lognormal) demand: average over what we don't know. The mixture is integrated by Gauss-Hermite quadrature (smooth tails). 'conf95' POS >= p when every rate sits at its 95th percentile. Reports how often that manifest really meets p: exactly 95% for a shared error, less when parts err separately. Convention: the datasheet rate is the lognormal's median, so the mean rate is exp(sigma^2 / 2) higher (about 9% at an error factor of 2). """ gh_z, gh_w = np.polynomial.hermite_e.hermegauss(80) gh_w = gh_w / gh_w.sum() rng = np.random.default_rng(seed) mu = demand_means(inp) mass = inp["mass"] nominal_cdfs = {i: (lambda k, m=mu[i]: stats.poisson.cdf(k, m)) for i in ITEMS} nominal = greedy_allocation(nominal_cdfs, mass, p) nominal_mass = sum(nominal[i] * mass[i] for i in ITEMS) out = {"design_target": p, "nominal_spares": nominal, "nominal_mass": nominal_mass, "by_ef": {}} for ef in efs: sig = lognormal_sigma(ef) if ef > 1 else 0.0 z_ind = rng.standard_normal((draws, len(ITEMS))) z_com = np.repeat(rng.standard_normal((draws, 1)), len(ITEMS), axis=1) res = {} for label, z in (("independent", z_ind), ("common", z_com)): x = np.exp(sig * z) pos = np.ones(draws) for j, i in enumerate(ITEMS): pos *= stats.poisson.cdf(nominal[i], mu[i] * x[:, j]) res[label] = { "pos_median": float(np.median(pos)), "pos_p05": float(np.quantile(pos, 0.05)), # 95% confidence level "pos_mean": float(pos.mean()), # = predictive POS "p_below_target": float((pos < p).mean()), } # hedged manifests if sig > 0: xs = np.exp(sig * gh_z) pred_cdfs = {i: (lambda k, m=mu[i], xs=xs: float(gh_w @ stats.poisson.cdf(k, m * xs))) for i in ITEMS} else: pred_cdfs = nominal_cdfs pred = greedy_allocation(pred_cdfs, mass, p) conf_cdfs = {i: (lambda k, m=mu[i] * ef: stats.poisson.cdf(k, m)) for i in ITEMS} conf = greedy_allocation(conf_cdfs, mass, p) res["predictive_spares"] = pred res["predictive_mass"] = sum(pred[i] * mass[i] for i in ITEMS) res["conf95_spares"] = conf res["conf95_mass"] = sum(conf[i] * mass[i] for i in ITEMS) for label, z in (("independent", z_ind), ("common", z_com)): x = np.exp(sig * z) pos = np.ones(draws) for j, i in enumerate(ITEMS): pos *= stats.poisson.cdf(conf[i], mu[i] * x[:, j]) res[f"conf95_meets_target_{label}"] = float((pos >= p).mean()) out["by_ef"][str(ef)] = res return out def aleatoric_vs_epistemic(inp=INPUTS, ef=2.0) -> dict: """Spread of demand for one item: Poisson scatter vs a rate known to EF.""" mu = demand_means(inp)["M1"] sig = lognormal_sigma(ef) pois_sd = math.sqrt(mu) # Var of Poisson-lognormal = mu E[X] + mu^2 Var[X] ex = math.exp(sig**2 / 2) varx = (math.exp(sig**2) - 1) * math.exp(sig**2) mixed_sd = math.sqrt(mu * ex + mu**2 * varx) return {"item": "M1", "mean": mu, "poisson_sd": pois_sd, "error_factor": ef, "mixed_sd": mixed_sd, "rate_part_of_variance": mu**2 * varx / (mu * ex + mu**2 * varx)} def learning(inp=INPUTS, failures=(3, 5, 10, 20, 30)) -> dict: """How well can testing pin down a rate? After r failures in exposure t, the Jeffreys posterior for the rate is Gamma(r + 1/2, t); its error factor is the 95th / 50th percentile ratio. Converted to years of the whole fleet running for M1 and the chassis at their nominal rates. """ n = inp["fleet"] rows = {} for r in failures: g = stats.gamma(r + 0.5) rows[str(r)] = {"error_factor": float(g.ppf(0.95) / g.ppf(0.5)), "fleet_years_M1": r * inp["core_mtbf_yr"] / inp["rate_share"]["M1"] / n, "fleet_years_chassis": r * inp["core_mtbf_yr"] / inp["rate_share"]["chassis"] / n} return rows # --- wear-out ----------------------------------------------------------------- def renewal_counts(n_slots, mean_life, beta, T, reps, rng, extra=None) -> np.ndarray: """Fleet failures in [0, T] for n_slots, each a renewal process. Lifetimes Weibull(beta) scaled to the given mean. If extra=(beta2, eta2) a second, independent wear-out mechanism competes (min of the two). Returns an array of length reps with the fleet total. """ eta = mean_life / math.gamma(1 + 1 / beta) totals = np.zeros(reps, dtype=np.int64) t = np.zeros((reps, n_slots)) alive = np.ones((reps, n_slots), dtype=bool) while alive.any(): life = eta * rng.weibull(beta, size=t.shape) if extra is not None: b2, e2 = extra life = np.minimum(life, e2 * rng.weibull(b2, size=t.shape)) t = np.where(alive, t + life, t) hit = alive & (t <= T) totals += hit.sum(axis=1) alive = hit return totals def wear_out(inp=INPUTS, p_item=0.9998, reps=50000, chunk=5000, seed=4077) -> dict: """Same MTBF, different ageing: what happens to M1's demand and spares? p_item is roughly the per-item level a 0.999 system target asks of the largest item (0.999 ** (1/5) ~ 0.9998). """ rng = np.random.default_rng(seed) n, T = inp["fleet"], inp["mission_yr"] mtbf_m1 = inp["core_mtbf_yr"] / inp["rate_share"]["M1"] cases = {} for label, beta, extra in ( ("infant (beta 0.7)", 0.7, None), ("random (beta 1)", 1.0, None), ("wear-out (beta 3)", 3.0, None), # random failures as specified, plus an unmodelled wear-out ("random + wear-out at eta 30 yr", 1.0, (4.0, 30.0)), ("random + wear-out at eta 45 yr", 1.0, (4.0, 45.0)), ): d = np.concatenate([renewal_counts(n, mtbf_m1, beta, T, chunk, rng, extra) for _ in range(reps // chunk)]) cases[label] = {"mean": float(d.mean()), "sd": float(d.std()), "spares": int(np.quantile(d, p_item, method="higher")), "poisson_spares_same_mean": poisson_quantile(float(d.mean()), p_item)} return {"item": "M1", "item_mtbf_yr": mtbf_m1, "p_item": p_item, "reps": reps, "cases": cases} # --- common cause ------------------------------------------------------------ def common_mode(inp=INPUTS, p=0.999, event_probs=(0.0, 0.001, 0.01, 0.05), hit_fracs=(0.1, 0.3)) -> dict: """A latent defect in one module type (a bad lot, a design flaw). With probability e over the mission an event fails a fraction q of the installed population of M1 at once, and the same fraction of the spares (they came from the same lot). Spares needed to ride it out: about q (n + s) extra. Report the POS of the nominal manifest with the event folded in, and what it would take to cover the event outright. """ mu = demand_means(inp) cdfs = {i: (lambda k, m=mu[i]: stats.poisson.cdf(k, m)) for i in ITEMS} s = greedy_allocation(cdfs, inp["mass"], p) base = pos_system(cdfs, s) n = inp["fleet"] others = math.prod(cdfs[i](s[i]) for i in ITEMS if i != "M1") rows = {} for q in hit_fracs: # given the event, M1 demand is D + Binomial(n + s, q); exact sum hits = stats.binom.pmf(np.arange(n + s["M1"] + 1), n + s["M1"], q) room = s["M1"] - np.arange(n + s["M1"] + 1) p_ok_event = float((hits * stats.poisson.cdf(room, mu["M1"])).sum()) * others # M1 spares to ride out the event at the per-item level (0.9998): # smallest s with P(D + Binomial(n + s, q) <= s) >= 0.9998 s_cover = s["M1"] while True: k = np.arange(n + s_cover + 1) ok = float((stats.binom.pmf(k, n + s_cover, q) * stats.poisson.cdf(s_cover - k, mu["M1"])).sum()) if ok >= 0.9998: break s_cover += 1 for e in event_probs: rows[f"q={q}, e={e}"] = {"pos_nominal_manifest": (1 - e) * base + e * p_ok_event, "pos_given_event": p_ok_event, "m1_spares_to_cover": s_cover, "extra_mass_cores": (s_cover - s["M1"]) * inp["mass"]["M1"], "m1_nominal": s["M1"]} return {"nominal_spares": s, "pos_without_event": base, "rows": rows} # --- figure data --------------------------------------------------------------- def demand_curves(inp=INPUTS, ef=2.0, nodes=200) -> dict: """M1 demand PMF with the rate known, and with the rate known to a factor EF. The mixed PMF averages Poisson(mu x) over a lognormal x by quadrature on equal-probability nodes. """ mu = demand_means(inp)["M1"] sig = lognormal_sigma(ef) x = np.exp(sig * stats.norm.ppf((np.arange(nodes) + 0.5) / nodes)) k = np.arange(0, 301) pois = stats.poisson.pmf(k, mu) mixed = stats.poisson.pmf(k[:, None], mu * x[None, :]).mean(axis=1) return {"k": k.tolist(), "poisson": pois.tolist(), "mixed": mixed.tolist(), "mean": mu, "error_factor": ef} def pos_vs_spares(inp=INPUTS) -> dict: """Fleet POS of manifest C as the M1 count varies, others at their 0.999 values.""" mu = demand_means(inp) cdfs = {i: (lambda k, m=mu[i]: stats.poisson.cdf(k, m)) for i in ITEMS} s = greedy_allocation(cdfs, inp["mass"], 0.999) others = math.prod(cdfs[i](s[i]) for i in ITEMS if i != "M1") ks = list(range(50, 121)) return {"m1": ks, "pos": [float(cdfs["M1"](k)) * others for k in ks], "m1_chosen": s["M1"]} def main() -> None: res = { "inputs": INPUTS, "strategies": strategies(), "k_of_n": k_of_n(), "aleatoric_vs_epistemic": aleatoric_vs_epistemic(), "rate_uncertainty": rate_uncertainty(), "learning": learning(), "wear_out": wear_out(), "common_mode": common_mode(), "demand_curves": demand_curves(), "pos_vs_spares": pos_vs_spares(), } (HERE / "results.json").write_text(json.dumps(res, indent=2)) print(json.dumps(res, indent=2)) if __name__ == "__main__": main()