MOSOF in ~100 lines of Python.

The MOSOF papers run on a few thousand lines of MATLAB and a Simulink-based fault simulator. The method underneath fits in about a hundred lines of Python. This is a walkthrough of those hundred lines — the inner loop, the fitness function, the Pareto-front extraction, and the knee selection — so you can read them in one sitting and know what is doing the work.

If you have read either of the Sensors (MDPI) papers,1Suslu, Ali, Jennions (2025); Suslu, Ali, Jennions (2026). Citations at the foot of this piece. you have seen the MOSOF block diagram: candidate pool in, multi-objective optimiser in the middle, Pareto front out, knee chosen by a stakeholder downstream. The block diagram earns the abstraction; the implementation has to make it actually run. The point of this walkthrough is to make the middle block legible — what the genetic algorithm is doing on each generation, what the three objectives evaluate to per suite, how the knee gets picked, and how all of that ends up in a 77-point Pareto front the rest of the system can render.

The full file lives in the mosof-ref-impl repository. The version below is the same code, broken into the sections that matter, with the rationale next to each.

01 The problem statement

MOSOF picks a subset of candidate sensors from a larger pool to optimise three competing objectives simultaneously:

• Diagnostic performance — NDCI-weighted coverage of the fault dictionary. Maximise.
• Suite cost — sum of unit costs across selected sensors. Minimise.
• Suite reliability — series-equivalent MTBF: 1 / Σ (1 / MTBF_i). Maximise.

The decision variable is a binary vector x ∈ {0,1}n where xi = 1 means sensor i is in the suite. The space is exponential in n — for a 32-sensor candidate pool, that is about 4 billion possible suites. Enumeration is out; we use a multi-objective genetic algorithm.

02 The candidate pool

The first thing you need is sensors to choose from. The reference impl seeds a pool whose per-family parameter ranges match the order-of-magnitude values in Section 4.3 of the thesis (the original per-sensor table is in the appendix; this is the synthetic-but-realistic version).

# Eight candidate sensors per subsystem, parameterised after
# thesis Section 4.3 family medians.
def make_candidate_pool(n_per_subsystem=8):
    rows = []
    for subsystem in ["Engine", "Fuel", "EPS", "ECS"]:
        cost_range = {"Engine": (3.0, 8.0), "Fuel": (1.5, 4.5),
                      "EPS": (2.0, 5.5), "ECS": (1.0, 3.5)}[subsystem]
        mtbf_range = {"Engine": (180, 480), "Fuel": (220, 520),
                      "EPS": (300, 700), "ECS": (240, 600)}[subsystem]
        for j in range(n_per_subsystem):
            rows.append({
                "id": f"{subsystem[:3].upper()}-{j+1:02d}",
                "subsystem": subsystem,
                "cost_kUSD": float(rng.uniform(*cost_range)),
                "mtbf_kh": float(rng.uniform(*mtbf_range)),
                "ndci": float(rng.uniform(0.05, 0.18)),
            })
    return rows

The four subsystems and their per-family cost / MTBF ranges aren't arbitrary — they follow the families used in the cross-subsystem study. The NDCI range [0.05, 0.18] is the per-sensor headline strength derived from the published fault × sensor matrix.2The /data/ndci-matrix.json dataset on the site exposes the published 12×18 matrix; the per-sensor headline strength is its column-sum, max-normalised.

03 The fitness function

For each candidate suite (a binary mask), compute the three objectives. This is the most opinionated part of the implementation — the choice of how coverage saturates and how reliability aggregates determines the shape of the front you get.

def evaluate_suite(mask, pool):
    sensors = pool["sensors"]
    idx = np.flatnonzero(mask.astype(bool))
    if idx.size == 0:
        return 0.0, 0.0, 0.0
    sel = [sensors[i] for i in idx]

    # Performance: NDCI-weighted coverage with diminishing returns.
    pool_ndci = sum(s["ndci"] for s in sensors)
    suite_ndci = sum(s["ndci"] for s in sel)
    perf = 1.0 - np.exp(-2.5 * suite_ndci / pool_ndci)

    # Cost: simple sum of unit costs.
    cost = sum(s["cost_kUSD"] for s in sel)

    # Reliability: series-equivalent MTBF.
    inv_mtbf = sum(1.0 / s["mtbf_kh"] for s in sel)
    mtbf = 1.0 / inv_mtbf if inv_mtbf > 0 else 0.0

    return float(perf), float(cost), float(mtbf)

Three choices to call out:

• The 1 - exp(-2.5 · ratio) form gives performance sub-additive saturation — the first sensors you add are worth a lot, the tenth is worth less, the twentieth almost nothing. This is what produces the knee. A linear coverage function would give you a front that is itself almost a straight line, and the rest of MOSOF would have nothing interesting to do.
• The cost is just a sum. Real installations have non-linear cost components (installation, certification, wiring run length), but those tend to be smooth functions of the count, and for a reference implementation the sum captures the essential trade-off without forcing the reader through a cost model they would have to learn to read.
• Series-equivalent MTBF — 1 / Σ(1 / MTBFi) — is the standard model for n components where any single failure breaks the diagnostic. Adding sensors strictly reduces suite MTBF; that is the third trade-off, and it is what stops "just add more sensors" from being a dominant strategy.

04 The NSGA-II driver

Now the actual optimisation. Pymoo's NSGA-II wants a problem class with an _evaluate method that returns the objective vector. We minimise, so the two maximise objectives get negated on the way out:

class MOSOFProblem(ElementwiseProblem):
    def __init__(self, pool):
        n = pool["n"]
        super().__init__(n_var=n, n_obj=3, n_ieq_constr=1,
                         xl=0, xu=1, vtype=int)
        self.pool = pool

    def _evaluate(self, x, out, *args, **kwargs):
        perf, cost, mtbf = evaluate_suite(x, self.pool)
        # NSGA-II minimises; negate the maximise axes.
        out["F"] = [-perf, cost, -mtbf]
        # Constraint: suite must have at least three sensors.
        out["G"] = [3 - int(x.sum())]

algorithm = NSGA2(
    pop_size=100,
    sampling=BinaryRandomSampling(),
    crossover=HalfUniformCrossover(),
    mutation=BitflipMutation(prob=1.0 / n),
    eliminate_duplicates=True,
)
result = minimize(MOSOFProblem(pool), algorithm, ("n_gen", 200), seed=1)

That is the whole optimisation. Three real choices worth knowing about:

• Half-uniform crossover over bit-flip mutation with rate 1/n. Standard NSGA-II for binary spaces. Half-uniform mixes bits that differ between parents; the 1/n mutation rate gives an expected one bit-flip per child, which is the long-time recommendation in the GA literature for this kind of selection problem.
• Population 100, 200 generations. Total budget is 20,000 fitness evaluations. For a 32-sensor pool this is more than enough; you can typically watch convergence stabilise around generation 80. The published runs use a larger budget; this is a reference impl, not a sweep.
• The constraint 3 - sum(x) ≤ 0 means "at least three sensors in the suite." NSGA-II treats feasible solutions as dominating infeasible ones, so the constraint never becomes a soft cost — it is a hard wall the front grows up against from below. The published MOSOF uses similar minimum-coverage constraints.

05 The Pareto front and the knee

Pymoo returns result.F already containing only non-dominated solutions — the algorithm's non-dominated sorting takes care of that. We do one more step: pick a knee.

def knee_by_perpendicular(front):
    perf = front[:, 0]
    cost = front[:, 1]
    perf_n = (perf - perf.min()) / max(perf.max() - perf.min(), 1e-9)
    cost_n = (cost - cost.min()) / max(cost.max() - cost.min(), 1e-9)
    p1 = np.array([perf_n[np.argmin(perf_n)],
                   cost_n[np.argmin(perf_n)]])
    p2 = np.array([perf_n[np.argmax(perf_n)],
                   cost_n[np.argmax(perf_n)]])
    line = p2 - p1
    line_len = np.linalg.norm(line) or 1e-9
    dists = []
    for i in range(len(front)):
        pt = np.array([perf_n[i], cost_n[i]])
        a, b = line, pt - p1
        d = (a[0] * b[1] - a[1] * b[0]) / line_len
        dists.append(abs(d))
    return int(np.argmax(dists))

This is the maximum-perpendicular-distance heuristic on the perf×cost projection: draw a line from the front's worst-performance endpoint to its best-performance endpoint, find the point farthest from that line. It is the simplest "elbow" detector and it works on most well-shaped fronts.

The published MOSOF uses a 3D measure that also weights reliability — pure 2D heuristics will pick a slightly different knee on the same front. That is a known property; the reference impl reports both, and downstream the canonical Pareto dataset on this site carries the published 3D knee explicitly in metadata so the page renders the right one.

The point of a reference implementation is to make the method legible — not to reproduce the published numbers down to the last decimal.

06 What you get back

Run it. On a laptop, NSGA-II converges in about five seconds. Output:

Candidate pool: 32 sensors across 4 subsystems

Pareto front: 100 non-dominated suites

Knee suite:
  size       = 17 sensors
  performance= 0.767  (normalised NDCI coverage)
  cost       = $50.5k
  MTBF       = 23 kh

Knee composition:
  Engine   2
  Fuel     3
  EPS      5
  ECS      7

Wrote pareto_front.csv (100 suites).

The published runs land closer to 12 sensors at performance ≈ 0.69 (Engine 5, Fuel 2, EPS 2, ECS 3) — that is the canonical knee.3Suslu, Ali, Jennions (2026), Section 4. The 12-sensor knee is the headline result of the cross-subsystem case. The reference impl's knee is larger because the synthetic pool's cost-to-NDCI ratio is more permissive than the published pool's. Tightening that ratio shrinks the knee back. The shape of the front is faithful; the exact knee composition depends on the pool you seed.

07 What this isn't

This is not the production MOSOF code. The production version has a high-fidelity Simulink-based fault simulator behind the NDCI evaluator, validation against thesis Table 4-5 with repeated nested cross-validation, sensor-failure models inside the reliability term, and per-subsystem constraints that the reference impl does not enforce. The papers describe what those add. The point of the reference impl is to make the method legible — what the GA is doing, what the three objectives are, how the knee gets picked — without the reader having to read the simulator first.

If you want to extend it: the most useful next step is replacing the synthetic NDCI per sensor with the real 12×18 matrix from /data/ndci-matrix.json and re-running. The ndci-calculator sibling repo already loads that matrix; lifting its ndci_per_sensor function into MOSOF's fitness gives you a front that lives much closer to the published numbers.

From there: the prior over fault subsystems becomes a real knob (the stakeholder-aware axis), the cost model can pull in per-family installation overheads, and the constraints can encode per-subsystem coverage minima. None of those changes the 100-line core. They change what the 100 lines are scoring.