Flow shop scheduling

Flow-Shop Scheduling

You are an expert in flow-shop scheduling. This skill covers the permutation flow-shop problem (PFSP): fast makespan evaluation, the NEH construction heuristic, exact MIP formulations, iterated greedy (IG) as the state-of-the-art heuristic, Taillard-style benchmark instances, and the multi-objective extension with makespan plus total tardiness. Use the framework below to match the formulation and solution method to instance size, objective, and time budget.

Initial Assessment

Establish the following before writing any model or code:

• Instance size. Number of jobs n and machines m. Taillard sizes run from 20x5 to 500x20. Exact methods prove optimality only for roughly n <= 15; everything larger is heuristic territory.
• Routing structure. Confirm every job visits machines 1..m in the same order. If jobs have job-specific routes, this is a job shop — switch to job-shop-scheduling.
• Permutation assumption. Is the same job order required on all machines? Standard PFSP assumes yes. For m >= 4, the best non-permutation schedule can beat the best permutation schedule; confirm the user accepts the (almost universal) permutation restriction.
• Buffer behavior. Unlimited intermediate buffers (standard), no buffers (blocking), or no waiting allowed between machines (no-wait)? Each changes the makespan recursion.
• Setup times. Sequence-dependent setup times (SDST) change both the MIP and the acceleration data structures; ask explicitly.
• Objective. Makespan C_max, total flowtime, total (weighted) tardiness, or several at once? Tardiness needs due dates — ask how they are generated or provided.
• Exact vs heuristic need. Is a provable optimum required, or is "within 1% of best known in 30 seconds" acceptable? This decides MIP vs NEH/IG immediately.
• Time budget. Per-instance wall-clock limit. The standard PFSP stopping convention is t = n * (m / 2) * 60 milliseconds (Ruiz & Stuetzle 2007); confirm what comparison protocol applies.
• Solver availability. Gurobi license present? If not, the MIP parts map to HiGHS/CBC via the same model structure, or skip exact methods entirely.
• Data format. Matrix of processing times p[i][j] (machine-major) or job-major? Taillard files list processing times machine by machine; off-by-one transposition is the most common data bug.
• Benchmark expectations. Will results be compared against Taillard or VRF best-known values? Then fix seeds, report relative percentage deviation (RPD), and pin the reference-value snapshot.
• Reproducibility. Seeds for instance generation and for every metaheuristic run; one result row per (instance, algorithm, seed).

Problem Definition and Model Landscape

Permutation flow-shop problem, Fm | prmu | C_max. Given n jobs and m machines, job j needs processing time p_{i,j} >= 0 on machine i. Every job visits machines 1, 2, ..., m in this order. All machines process jobs in the same sequence pi (a permutation of the jobs). Each machine handles one job at a time; operations are non-preemptive. Completion times follow the recursion

$$ C_{i,\pi(k)} = \max\big(C_{i-1,\pi(k)},; C_{i,\pi(k-1)}\big) + p_{i,\pi(k)}, \qquad C_{0,\cdot} = 0,\quad C_{\cdot,\pi(0)} = 0, $$

for machines i = 1..m and sequence positions k = 1..n. The makespan is C_max(pi) = C_{m, pi(n)}. The search space is the set of n! permutations.

Complexity landmarks:

• F2 || C_max is solvable in O(n log n) by Johnson's rule (Johnson 1954, "Optimal two- and three-stage production schedules with setup times included").
• F3 || C_max is already strongly NP-hard (Garey, Johnson & Sethi 1976, "The complexity of flowshop and jobshop scheduling").
• For m <= 3 an optimal permutation schedule exists; for m >= 4 the gap between the best permutation and the best non-permutation schedule can grow with m (Potts, Shmoys & Williamson 1991). In practice almost all research and industry use restricts to permutations.

Variants

Variant	What changes	Practical note
Standard PFSP Fm|prmu|C_max	nothing — the base case	this skill's core
No-wait Fm|nwt|C_max	a job may never wait between machines	reduces to an asymmetric TSP over completion-distance
Blocking Fm|block|C_max	no intermediate buffers; a finished job blocks its machine	different recursion (departure times), different acceleration
No-idle	machines must run without idle time	niche; changes feasibility of sequences
SDST	sequence-dependent setup s_{i,j,k} between job j and k on machine i	add setups to recursion and MIP; IG still works well
Hybrid / flexible flow shop	parallel machines per stage	assignment + sequencing; combine with parallel-machine ideas
Distributed PFSP	several factories, assign then sequence	two-level decisions; IG variants dominate

Objectives

Objective	Definition	Method of choice
Makespan C_max	completion of last job on last machine	NEH + iterated greedy
Total flowtime sum C_j	sum of last-machine completions	IG variants with different acceleration (Pan & Ruiz 2012)
Total tardiness sum T_j	T_j = max(C_j - d_j, 0) with due dates d_j	IG/local search; due-date generation matters
Bi-objective (C_max, sum T_j)	Pareto front	weighted sum / epsilon-constraint / NSGA-II — see multi-objective-optimization

Method selection

• n <= 12-15: solve the MIP to optimality; warm start with NEH. Use the exact optimum to validate the heuristic pipeline.
• Any n, milliseconds available: NEH. It lands within roughly 2-4% of best known on Taillard instances and is the standard starting point for everything else.
• Seconds to minutes available: iterated greedy. Despite hundreds of published "novel" metaheuristics, IG and close relatives remain the state of the art for PFSP makespan (Ruiz & Stuetzle 2007; Fernandez-Viagas, Ruiz & Framinan 2017 review).
• Population-based alternative (e.g., when a decoder-based framework already exists): BRKGA with a sort decoder — see biased-random-key-genetic-algorithm.
• Non-permutation, blocking with complex side constraints, or mixed shop features: consider the CP-SAT interval model from job-shop-scheduling instead of forcing a PFSP shape.

Benchmarks

Taillard (1993, "Benchmarks for basic scheduling problems"): 120 instances, 12 sizes (20x5 up to 500x20), 10 instances per size, processing times integer uniform on [1, 99] from a published seed list. Vallada, Ruiz & Framinan (2015) add 480 larger, harder instances (VRF). Report performance as relative percentage deviation from the best-known value: RPD = 100 * (C - C_best) / C_best, averaged over instances and seeds.

Makespan Evaluation and the NEH Heuristic

Everything in PFSP work stands on a fast, correct makespan function. The recursion is a max-plus scan that is sequential in both dimensions, so the right vectorization axis is the batch of permutations (populations, multi-seed evaluation), not a single sequence.

import numpy as np


def makespan(p: np.ndarray, perm: np.ndarray) -> int:
    """C_max of permutation `perm` for processing times p with shape (m, n).

    Rolling-array form of the recursion; O(nm) time, O(n) memory. The max-plus
    scan is sequential in both dimensions, so these two loops cannot be replaced
    by broadcasting; vectorize over batches of permutations instead.
    """
    q = p[:, perm]
    c = np.zeros(q.shape[1], dtype=np.int64)
    for i in range(q.shape[0]):
        c[0] += q[i, 0]
        for k in range(1, q.shape[1]):
            c[k] = max(c[k], c[k - 1]) + q[i, k]
    return int(c[-1])


def completion_matrix(p: np.ndarray, perm: np.ndarray) -> np.ndarray:
    """Full (m, n) completion-time matrix in sequence order.

    C[i, k] is the completion of the k-th sequenced job on machine i; use it
    for Gantt charts, tardiness objectives, and validators."""
    m, n = p.shape
    q = p[:, perm]
    c = np.zeros((m, n), dtype=np.int64)
    for i in range(m):
        for k in range(n):
            left = c[i, k - 1] if k > 0 else 0
            up = c[i - 1, k] if i > 0 else 0
            c[i, k] = max(left, up) + q[i, k]
    return c


def makespan_batch(p: np.ndarray, perms: np.ndarray) -> np.ndarray:
    """Makespans of a (B, n) batch of permutations, vectorized over the batch.

    Use this to evaluate whole metaheuristic populations in one call."""
    m = p.shape[0]
    q = p[:, perms]                       # (m, B, n) via fancy indexing
    c = np.zeros(perms.shape, dtype=np.int64)
    for i in range(m):
        c[:, 0] += q[i, :, 0]
        for k in range(1, perms.shape[1]):
            c[:, k] = np.maximum(c[:, k], c[:, k - 1]) + q[i, :, k]
    return c[:, -1]


p = np.array([[3, 2, 4], [2, 4, 1]])      # m=2 machines, n=3 jobs
print(makespan(p, np.array([0, 1, 2])))   # Expected: 10 (optimal; Johnson's rule confirms)
print(makespan_batch(p, np.array([[0, 1, 2], [2, 0, 1]])))
# Expected: [10 13]

NEH (Nawaz, Enscore & Ham 1983) sorts jobs by decreasing total processing time and inserts them one at a time at the makespan-minimizing position. Implemented naively it costs O(n^3 m); with Taillard's head/tail acceleration (Taillard 1990, "Some efficient heuristic methods for the flow shop sequencing problem") all k+1 insertion positions of one job are evaluated in O(mk) total, giving O(n^2 m) for the whole heuristic. The same routine is the engine inside iterated greedy.

For a partial sequence s of length k, define heads e[i, l] (earliest completion of s[l] on machine i), tails t[i, l] (time from the start of s[l] on machine i to the end of the schedule, computed backwards), and f[i, l] (completion of the inserted job at position l on machine i). The makespan after inserting at position l is max_i (f[i, l] + t[i, l]), with a zero tail when inserting at the end.

import numpy as np


def insertion_makespans(p: np.ndarray, seq: np.ndarray, job: int) -> np.ndarray:
    """Makespan of inserting `job` at every position 0..len(seq) of `seq`.

    Taillard (1990) acceleration: heads e, tails t, and inserted-job completions
    f give all len(seq)+1 makespans in O(m * len(seq)) total, instead of
    O(m * len(seq)^2) for naive re-evaluation of every position.
    """
    m, k = p.shape[0], seq.size
    q = p[:, seq].astype(np.float64)
    e = np.zeros((m, k))       # e[i, l]: completion of seq[l] on machine i
    t = np.zeros((m, k))       # t[i, l]: tail from the start of seq[l] on machine i
    f = np.zeros((m, k + 1))   # f[i, l]: completion of `job` inserted at position l
    for i in range(m):
        for l in range(k):
            e[i, l] = max(e[i, l - 1] if l else 0.0,
                          e[i - 1, l] if i else 0.0) + q[i, l]
    for i in range(m - 1, -1, -1):
        for l in range(k - 1, -1, -1):
            t[i, l] = max(t[i, l + 1] if l < k - 1 else 0.0,
                          t[i + 1, l] if i < m - 1 else 0.0) + q[i, l]
    for l in range(k + 1):
        for i in range(m):
            f[i, l] = max(e[i, l - 1] if l else 0.0,
                          f[i - 1, l] if i else 0.0) + p[i, job]
    tails = np.hstack([t, np.zeros((m, 1))])   # inserting at the end has no tail
    return (f + tails).max(axis=0)


def neh(p: np.ndarray) -> tuple[np.ndarray, int]:
    """NEH heuristic (Nawaz, Enscore & Ham 1983), O(n^2 m) with acceleration.

    Sort jobs by decreasing total processing time (stable, so ties keep index
    order — document this, ties change the result), then insert each job at the
    position minimizing the partial-sequence makespan.
    """
    order = np.argsort(-p.sum(axis=0), kind="stable")
    seq = order[:1].copy()
    best = int(p[:, order[0]].sum())
    for job in order[1:]:
        ms = insertion_makespans(p, seq, int(job))
        pos = int(np.argmin(ms))
        best = int(ms[pos])
        seq = np.insert(seq, pos, job)
    return seq, best


p = np.array([[3, 2, 4], [2, 4, 1]])
perm, cmax = neh(p)
print(perm, cmax)
# Expected: perm [1 0 2] with makespan 10 (NEH finds the optimum on this instance)

Exact MIP Models in Gurobi

Two classic formulation families exist (computational comparison: Tseng, Stafford & Gupta 2004). Both prove optimality only on small instances; their real value is validating heuristics and producing certified optima/bounds for papers.

Positional model (Wilson 1989 family). Binary x[j,k] = 1 if job j occupies sequence position k; continuous c[i,k] is the completion time of the position-k job on machine i:

$$ \min\ c_{m,n} \quad \text{s.t.}\quad \sum_k x_{jk} = 1,\ \sum_j x_{jk} = 1;\qquad c_{ik} \ge c_{i-1,k} + \textstyle\sum_j p_{ij} x_{jk};\qquad c_{ik} \ge c_{i,k-1} + \textstyle\sum_j p_{ij} x_{jk}. $$

No big-M appears: position-indexed completion times encode both the machine order and the sequence order directly, which gives a comparatively strong LP relaxation at the price of n^2 binaries.

import gurobipy as gp
import numpy as np
from gurobipy import GRB


def add_assignment_constraints(model: gp.Model, vars_: dict, data: np.ndarray) -> None:
    """Each job takes exactly one position; each position holds exactly one job."""
    x = vars_["x"]
    n = data.shape[1]
    model.addConstrs((x.sum(j, "*") == 1 for j in range(n)), name="job_once")
    model.addConstrs((x.sum("*", k) == 1 for k in range(n)), name="pos_once")


def add_route_constraints(model: gp.Model, vars_: dict, data: np.ndarray) -> None:
    """The job in position k starts on machine i only after leaving machine i-1."""
    x, c = vars_["x"], vars_["c"]
    m, n = data.shape
    for i in range(m):
        for k in range(n):
            proc = gp.quicksum(float(data[i, j]) * x[j, k] for j in range(n))
            prev = c[i - 1, k] if i > 0 else 0.0
            model.addConstr(c[i, k] >= prev + proc, name=f"route[{i},{k}]")


def add_sequence_constraints(model: gp.Model, vars_: dict, data: np.ndarray) -> None:
    """Machine i starts position k only after finishing position k-1."""
    x, c = vars_["x"], vars_["c"]
    m, n = data.shape
    for i in range(m):
        for k in range(1, n):
            proc = gp.quicksum(float(data[i, j]) * x[j, k] for j in range(n))
            model.addConstr(c[i, k] >= c[i, k - 1] + proc, name=f"seq[{i},{k}]")


def solve_pfsp_positional(p: np.ndarray, time_limit: float = 60.0,
                          warm_start: np.ndarray | None = None) -> tuple[list[int], float, float]:
    """Positional PFSP MIP. Returns (permutation, C_max, final MIP gap)."""
    m, n = p.shape
    model = gp.Model("pfsp_positional")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    x = model.addVars(n, n, vtype=GRB.BINARY, name="x")
    c = model.addVars(m, n, lb=0.0, name="c")
    vars_ = {"x": x, "c": c}
    add_assignment_constraints(model, vars_, p)
    add_route_constraints(model, vars_, p)
    add_sequence_constraints(model, vars_, p)
    model.setObjective(c[m - 1, n - 1], GRB.MINIMIZE)
    if warm_start is not None:               # e.g. the NEH permutation (MIP start)
        for k, j in enumerate(warm_start):
            x[int(j), int(k)].Start = 1.0
    model.optimize()
    ok = model.Status == GRB.OPTIMAL or (model.Status == GRB.TIME_LIMIT and model.SolCount > 0)
    if not ok:
        raise RuntimeError(f"no feasible solution, status {model.Status}")
    perm = [j for k in range(n) for j in range(n) if x[j, k].X > 0.5]
    return perm, model.ObjVal, model.MIPGap


p = np.array([[3, 2, 4], [2, 4, 1]])
perm, cmax, gap = solve_pfsp_positional(p, time_limit=10.0)
print(perm, cmax, gap)
# Expected: a permutation with C_max = 10.0 and gap 0.0 (e.g. [0, 1, 2])

Disjunctive model (Manne 1960 family). One precedence binary y[j,l] per unordered job pair, shared across all machines — sharing is exactly what enforces the permutation property. Big-M constraints order each pair on each machine. Only n(n-1)/2 binaries, but the LP relaxation is weak because of the big-M terms. In the Tseng-Stafford-Gupta experiments Manne-type models are often the fastest to solve despite the weaker bound, because they are much smaller — test both on your instances.

import gurobipy as gp
import numpy as np
from gurobipy import GRB


def add_route_constraints(model: gp.Model, vars_: dict, data: np.ndarray) -> None:
    """Each job follows the machine order 1..m."""
    c = vars_["c"]
    m, n = data.shape
    for j in range(n):
        model.addConstr(c[0, j] >= float(data[0, j]), name=f"first[{j}]")
        for i in range(1, m):
            model.addConstr(c[i, j] >= c[i - 1, j] + float(data[i, j]),
                            name=f"route[{i},{j}]")


def add_disjunctive_constraints(model: gp.Model, vars_: dict, data: np.ndarray) -> None:
    """One precedence binary per job pair, shared by all machines (permutation property)."""
    c, y = vars_["c"], vars_["y"]
    m, n = data.shape
    big_m = float(data.sum())            # tightest simple bound: total work
    for j in range(n):
        for l in range(j + 1, n):
            for i in range(m):
                model.addConstr(c[i, j] >= c[i, l] + float(data[i, j]) - big_m * y[j, l],
                                name=f"disj_a[{i},{j},{l}]")
                model.addConstr(c[i, l] >= c[i, j] + float(data[i, l]) - big_m * (1 - y[j, l]),
                                name=f"disj_b[{i},{j},{l}]")


def add_makespan_constraints(model: gp.Model, vars_: dict, data: np.ndarray) -> None:
    """C_max dominates every job's completion on the last machine."""
    c, cmax = vars_["c"], vars_["cmax"]
    m, n = data.shape
    model.addConstrs((cmax >= c[m - 1, j] for j in range(n)), name="span")


def solve_pfsp_manne(p: np.ndarray, time_limit: float = 60.0) -> tuple[list[int], float, float]:
    """Manne-style disjunctive PFSP model: n(n-1)/2 binaries, big-M precedence."""
    m, n = p.shape
    model = gp.Model("pfsp_manne")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    c = model.addVars(m, n, lb=0.0, name="c")
    y = model.addVars(((j, l) for j in range(n) for l in range(j + 1, n)),
                      vtype=GRB.BINARY, name="y")
    cmax = model.addVar(lb=0.0, name="cmax")
    vars_ = {"c": c, "y": y, "cmax": cmax}
    add_route_constraints(model, vars_, p)
    add_disjunctive_constraints(model, vars_, p)
    add_makespan_constraints(model, vars_, p)
    model.setObjective(cmax, GRB.MINIMIZE)
    model.optimize()
    ok = model.Status == GRB.OPTIMAL or (model.Status == GRB.TIME_LIMIT and model.SolCount > 0)
    if not ok:
        raise RuntimeError(f"no feasible solution, status {model.Status}")
    perm = sorted(range(n), key=lambda j: c[m - 1, j].X)
    return perm, cmax.X, model.MIPGap


p = np.array([[3, 2, 4], [2, 4, 1]])
perm, cmax_val, gap = solve_pfsp_manne(p, time_limit=10.0)
print(perm, cmax_val, gap)
# Expected: C_max = 10.0 with gap 0.0; perm read off last-machine completion order

Feed the NEH permutation into warm_start whenever you run the positional model — the solver gets an incumbent immediately and the search becomes pure bound-proving. Patterns for MIP starts and variable hints live in warm-starts-and-initial-solutions.

Iterated Greedy — the State of the Art

Iterated greedy (Ruiz & Stuetzle 2007, "A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem") repeats two steps from an NEH start: destruction (remove d random jobs) and construction (reinsert each removed job at its best position with the Taillard acceleration), with a constant-temperature Metropolis acceptance. It is the IG instantiation of the ILS template — perturbation-strength tuning, acceptance variants, and restart logic are covered in iterated-local-search. The insertion evaluation is the vectorizable hot spot: each call scores all k+1 positions in one accelerated pass.

import numpy as np


def _makespan(p: np.ndarray, perm: np.ndarray) -> float:
    """C_max via the rolling-array recursion; O(nm)."""
    q = p[:, perm]
    c = np.zeros(q.shape[1])
    for i in range(q.shape[0]):
        c[0] += q[i, 0]
        for k in range(1, q.shape[1]):
            c[k] = max(c[k], c[k - 1]) + q[i, k]
    return float(c[-1])


def _best_insertion(p: np.ndarray, seq: np.ndarray, job: int) -> tuple[int, float]:
    """Best insertion position and resulting C_max (Taillard acceleration, O(m*len(seq)))."""
    m, k = p.shape[0], seq.size
    q = p[:, seq].astype(np.float64)
    e = np.zeros((m, k))
    t = np.zeros((m, k))
    f = np.zeros((m, k + 1))
    for i in range(m):
        for l in range(k):
            e[i, l] = max(e[i, l - 1] if l else 0.0, e[i - 1, l] if i else 0.0) + q[i, l]
    for i in range(m - 1, -1, -1):
        for l in range(k - 1, -1, -1):
            t[i, l] = max(t[i, l + 1] if l < k - 1 else 0.0,
                          t[i + 1, l] if i < m - 1 else 0.0) + q[i, l]
    for l in range(k + 1):
        for i in range(m):
            f[i, l] = max(e[i, l - 1] if l else 0.0, f[i - 1, l] if i else 0.0) + p[i, job]
    ms = (f + np.hstack([t, np.zeros((m, 1))])).max(axis=0)
    pos = int(np.argmin(ms))
    return pos, float(ms[pos])


def iterated_greedy(p: np.ndarray, seed: int = 0, max_iters: int = 2000, d: int = 4,
                    temp_factor: float = 0.4) -> tuple[np.ndarray, float]:
    """IG of Ruiz & Stuetzle (2007): NEH start, destruct-reconstruct loop,
    constant-temperature acceptance. Framework: see iterated-local-search."""
    rng = np.random.default_rng(seed)
    m, n = p.shape
    temp = temp_factor * p.sum() / (10.0 * n * m)
    order = np.argsort(-p.sum(axis=0), kind="stable")    # NEH construction
    cur = order[:1].copy()
    for job in order[1:]:
        pos, _ = _best_insertion(p, cur, int(job))
        cur = np.insert(cur, pos, job)
    cur_ms = _makespan(p, cur)
    best, best_ms = cur.copy(), cur_ms
    for _ in range(max_iters):
        removed = rng.choice(n, size=min(d, n - 1), replace=False)
        partial = cur[~np.isin(cur, removed)]
        cand_ms = cur_ms
        for job in removed:                              # greedy reconstruction
            pos, cand_ms = _best_insertion(p, partial, int(job))
            partial = np.insert(partial, pos, job)
        if cand_ms < cur_ms or rng.random() < np.exp((cur_ms - cand_ms) / temp):
            cur, cur_ms = partial, cand_ms
            if cur_ms < best_ms:
                best, best_ms = cur.copy(), cur_ms
    return best, best_ms


p = np.array([[3, 2, 4], [2, 4, 1]], dtype=np.float64)
perm, cmax = iterated_greedy(p, seed=42, max_iters=200, d=2)
print(perm, cmax)
# Expected: C_max = 10.0 (the optimum for this 3-job, 2-machine instance)

Parameter guidance (defaults are robust — IG is famously insensitive):

Parameter	Typical range	Trade-off
d (jobs removed)	2-8, default 4	larger = stronger perturbation, slower reconstruction
temp_factor	0.2-0.7, default 0.4	higher accepts more worsening moves (diversification)
stopping	n*(m/2)*60 ms or fixed iterations	time-based for benchmark comparability

The full IG_RS adds an insertion local search (repeatedly extract one job, reinsert at its best position, until a full pass yields no improvement) after each reconstruction; it costs O(n^2 m) per pass and buys roughly 0.3-0.6% average RPD on Taillard instances. Add it when the time budget allows more than a few seconds.

Instance Generation and Independent Validation

Group every experiment around seeded generators and a validator that shares no code with the optimization pipeline. The generator below follows Taillard's distribution; the due-date scheme follows the tardiness literature (Vallada, Ruiz & Minella 2008, tardiness benchmark design).

import numpy as np


def taillard_instance(m: int, n: int, seed: int) -> np.ndarray:
    """Taillard (1993)-style PFSP instance: integer processing times U[1, 99].

    Returns p with shape (m, n). Fix the seed and report (m, n, seed) next to
    every result row so any number in a paper can be regenerated."""
    rng = np.random.default_rng(seed)
    return rng.integers(1, 100, size=(m, n), dtype=np.int64)


def lower_bound(p: np.ndarray) -> int:
    """Machine-based lower bound on C_max (Taillard 1993).

    For each machine: its total load plus the smallest head (work before it)
    and smallest tail (work after it) over jobs; also the longest job."""
    m, n = p.shape
    heads = np.vstack([np.zeros(n), np.cumsum(p, axis=0)[:-1]])   # work before machine i
    tails = p[::-1].cumsum(axis=0)[::-1] - p                      # work after machine i
    machine_lb = (p.sum(axis=1) + heads.min(axis=1) + tails.min(axis=1)).max()
    job_lb = p.sum(axis=0).max()
    return int(max(machine_lb, job_lb))


def due_dates(p: np.ndarray, seed: int, tf: float = 0.4, rdd: float = 0.6) -> np.ndarray:
    """Due dates for tardiness objectives via the TF/RDD scheme.

    Drawn uniformly from [LB*(1-tf-rdd/2), LB*(1-tf+rdd/2)] where LB estimates
    C_max. Larger tf = tighter dates = more tardy jobs; larger rdd = wider
    spread. Both shift instance difficulty, so always report (tf, rdd, seed)."""
    rng = np.random.default_rng(seed)
    lb = lower_bound(p)
    lo, hi = lb * (1.0 - tf - rdd / 2.0), lb * (1.0 - tf + rdd / 2.0)
    d = rng.uniform(max(lo, 0.0), hi, size=p.shape[1])
    return np.maximum(d, p.sum(axis=0))    # a job is never due before its own total work


p_small = np.array([[3, 2, 4], [2, 4, 1]])
print(lower_bound(p_small))    # Expected: 10 — matches the optimal C_max on this instance
p = taillard_instance(5, 20, seed=1)
print(p.shape, bool(p.min() >= 1 and p.max() <= 99))
# Expected: (5, 20) True

The validator recomputes feasibility and the objective in plain Python, on purpose: a bug in the vectorized evaluator (or in the MIP) cannot also live in this code path. Run it on every solution that leaves the pipeline, and cross-check exact-vs-heuristic results on small instances.

import numpy as np


def validate_pfsp(p: np.ndarray, perm: np.ndarray, claimed_cmax: float,
                  tol: float = 1e-6) -> list[str]:
    """Independent feasibility + objective check; returns violations (empty == valid).

    Checks: (1) perm is a true permutation of 0..n-1; (2) the claimed makespan
    equals a from-scratch completion-time recomputation done in plain Python,
    independent of any vectorized or solver code."""
    m, n = p.shape
    if sorted(int(j) for j in perm) != list(range(n)):
        return [f"perm is not a permutation of 0..{n - 1}: {list(perm)}"]
    c = [[0.0] * n for _ in range(m)]
    for k, j in enumerate(int(v) for v in perm):
        for i in range(m):
            left = c[i][k - 1] if k > 0 else 0.0
            up = c[i - 1][k] if i > 0 else 0.0
            c[i][k] = max(left, up) + float(p[i, j])
    recomputed = c[m - 1][n - 1]
    if abs(recomputed - float(claimed_cmax)) > tol:
        return [f"claimed C_max {claimed_cmax} != recomputed {recomputed}"]
    return []


p = np.array([[3, 2, 4], [2, 4, 1]])
print(validate_pfsp(p, np.array([1, 0, 2]), 10))
# Expected: [] (valid)
print(validate_pfsp(p, np.array([1, 0, 2]), 9))
# Expected: ['claimed C_max 9 != recomputed 10.0']

Advanced Techniques

Taillard acceleration beyond NEH

The head/tail data structure pays off anywhere insertions are evaluated: NEH (O(n^3 m) down to O(n^2 m)), IG reconstruction, and the insertion local search. Two caveats. First, the trick is makespan-specific — tails do not exist for total flowtime or tardiness, where every completion time matters; flowtime needs different incremental schemes (Pan & Ruiz 2012) or honest recomputation. Second, heads e can be reused across the d insertions of one IG reconstruction only for the prefix untouched by earlier insertions; the simple version above recomputes them, which is still the standard implementation and rarely the bottleneck below n = 500.

NEH tie-breaking and priority orders

Two tie sources change NEH results between implementations: equal job totals in the initial sort and equal makespans in argmin over positions. A stable sort plus first-best argmin (as coded above) is reproducible but not best-performing. NEH-D (Dong, Huang & Chen 2008) orders jobs by mean plus standard deviation of processing times; the tie-breaker of Fernandez-Viagas & Framinan (2014) picks among equal-makespan positions the one minimizing front-machine idle time and is the strongest published variant. When you compare against published NEH numbers, state the exact tie-breaking rule — differences up to 1% RPD come from ties alone.

Tuning and extending iterated greedy

IG's strength is robustness: on Taillard benchmarks, d = 4 and temp_factor = 0.4 are near-optimal across all sizes, so tune the stopping rule before anything else. Worthwhile extensions, in order of payoff: (1) insertion local search after reconstruction (IG_RS); (2) local search also on the partial sequence after destruction (IG variants of Fernandez-Viagas & Framinan); (3) for total flowtime or tardiness, swap the acceptance to best-insertion under the new objective and replace the acceleration. Population-based methods (BRKGA with a sort decoder, EDAs with position-frequency models) reach similar quality but need more evaluations; they pay off when a decoder framework or parallel fitness evaluation already exists.

Bi-objective makespan + tardiness

Bi-objective PFSP (C_max, sum T_j) is the standard multi-objective extension (survey: Minella, Ruiz & Ciavotta 2008). For a handful of trade-off points run epsilon-constraint on the MIP (constrain sum T_j <= eps, minimize C_max, sweep eps); for a dense front use NSGA-II or a Pareto-archive IG. Non-dominated sorting, crowding, hypervolume, and pymoo usage are in multi-objective-optimization. The building blocks:

import numpy as np


def bi_objective(p: np.ndarray, d: np.ndarray, perm: np.ndarray) -> tuple[float, float]:
    """(C_max, total tardiness) of a permutation; d[j] is job j's due date."""
    q = p[:, perm]
    c = np.zeros(q.shape[1])
    for i in range(q.shape[0]):
        c[0] += q[i, 0]
        for k in range(1, q.shape[1]):
            c[k] = max(c[k], c[k - 1]) + q[i, k]
    tardiness = float(np.maximum(c - d[perm], 0.0).sum())
    return float(c[-1]), tardiness


def pareto_front(points: np.ndarray) -> np.ndarray:
    """Indices of non-dominated rows of `points` (minimization), vectorized.

    Duplicates are kept: identical points do not dominate each other."""
    le = (points[:, None, :] <= points[None, :, :]).all(axis=2)
    lt = (points[:, None, :] < points[None, :, :]).any(axis=2)
    dominated = (le & lt).any(axis=0)
    return np.flatnonzero(~dominated)


p = np.array([[3, 2, 4], [2, 4, 1]])
d = np.array([6.0, 5.0, 11.0])
pts = np.array([bi_objective(p, d, np.array(s))
                for s in ([0, 1, 2], [1, 0, 2], [2, 1, 0])])
print(pts[pareto_front(pts)])
# Expected: [[10. 3.]] — perm [1, 0, 2] dominates the other two sequences

Lower bounds and gap reporting

The machine-based bound in lower_bound is cheap and surprisingly tight on Taillard instances with n >> m (often within 1-2% of the optimum for 100x5). Use it three ways: sanity-check every reported makespan (C >= LB always), report heuristic gaps on instances without best-known values (100 * (C - LB) / LB is then an upper bound on the true gap), and prune in any custom branch-and-bound. The LP relaxation of the positional model gives a stronger but far more expensive bound; it is rarely worth solving for bounds alone.

Practical Challenges

The MIP stalls beyond about 15 jobs. This is expected, not a modeling bug: PFSP MIPs have weak relaxations relative to the combinatorial structure, and n = 20 is already out of reach for proof of optimality in reasonable time. Use the MIP to certify optima on small instances (validating the heuristic pipeline), warm start it with NEH, and report the residual MIPGap honestly when hitting the time limit.

Your IG results are worse than published numbers. Check, in order: time-based stopping with the n*(m/2)*60 ms convention versus your iteration cap; presence of the insertion local search (published IG_RS includes it); the temperature constant (the denominator 10*n*m is part of the formula, not a typo); and the RPD protocol (averaging over the right instance set, 5+ seeds per instance, best-known values from the same snapshot the paper used).

NEH gives different makespans across implementations. Ties in the initial sort and in the insertion argmin are resolved differently. Fix a deterministic rule (stable sort, first minimizing position), document it, and never compare NEH variants across codebases without matching tie-breakers.

Makespan evaluation dominates the runtime. Profile before optimizing, but the usual fixes are: never rebuild the full completion matrix when the accelerated insertion routine answers the same question; evaluate populations with makespan_batch instead of a Python loop over individuals; and keep processing times in one contiguous (m, n) int64 array — repeated fancy indexing of ragged structures costs more than the recursion itself.

The permutation assumption silently discards better schedules. For m >= 4 the best non-permutation schedule can beat the best permutation one. Almost all benchmarks and applications accept this restriction, but say so explicitly in any write-up; if buffers are actually finite or zero, the model must change (blocking PFSP), not just the search space.

Tardiness results are not comparable across papers. Due-date generation parameters (TF, RDD) move instances between trivially easy (everything early) and brutally tight. Publish the generator, the parameters, and the seeds with the results, and never mix due-date schemes inside one experiment table.

Best-known values drift over the years. Taillard upper bounds have been improved repeatedly since 1993. Pin the snapshot of best-known values you compare against (date + source), or your RPD numbers cannot be reproduced after the next improvement.

Integer data, float drift. Taillard processing times are integers, so makespans are integers; compute and compare them in int64 where possible. The float versions inside IG are fine for search, but validate final results with exact integer recomputation to avoid reporting a 1245.0000000002.

Tools & Libraries

Library	When to use	Note
gurobipy	exact MIP models, warm-started optimality proofs	both formulations above; commercial license
OR-Tools CP-SAT	shop variants with side constraints, non-permutation schedules	interval variables; see job-shop-scheduling
HiGHS (highspy) / PuLP	no Gurobi license	same model structure; expect smaller solvable sizes
numpy	makespan evaluation, NEH, IG, batch population evaluation	the whole heuristic stack needs nothing else
pandas	result tables: one row per (instance, algorithm, seed)	aggregation to mean/best RPD per size
matplotlib	Gantt charts from completion_matrix, convergence curves	vector output for papers
pymoo	NSGA-II and quality indicators for the bi-objective extension	pairs with the bi_objective evaluator

Output Format

A complete flow-shop answer contains:

1. Problem statement echo. n, m, objective(s), variant (standard/blocking/SDST), permutation assumption confirmed, data source (Taillard set or generator with seeds).

2. Model/method summary table.

   Component	Choice	Why
   Construction	NEH (stable sort, first-best insertion)	standard, reproducible baseline
   Improvement	IG, d=4, temp_factor=0.4, time n*(m/2)*60 ms	state of the art for C_max
   Exact reference	positional MIP, NEH warm start, 300 s limit	optima for n <= 12 validation set
   Validation	validate_pfsp on every reported solution	independent recomputation

3. Solution-quality report. Per instance size: best/mean/worst RPD over seeds, mean time to best, and for MIP runs the bound, incumbent, and final gap. State the best-known-value snapshot.

4. Validator confirmation. Explicit statement that every reported permutation passed the independent check, plus the exact-vs-heuristic cross-check result on the small-instance set.

5. Artifacts. results.csv (one row per run: instance, m, n, seed, algorithm, parameters, objective, time), instance files or generator seeds, and optionally gantt.png of the best schedule and a convergence plot (best-so-far C_max vs time, band over seeds).

6. Reproduction line. The exact command/function call and seed list that regenerate every number shown.

Questions to Ask

• How many jobs and machines, and is the size fixed or growing in production use?
• Does every job really visit all machines in the same order, or do some skip machines or take different routes?
• Is the same sequence on all machines acceptable (permutation schedules), or do you need the extra freedom of non-permutation schedules?
• Are there buffers between machines — unlimited, limited, or none (blocking)?
• Are setup times relevant, and do they depend on the job sequence?
• Which objective: makespan, flowtime, tardiness, or a trade-off between several?
• If tardiness matters: where do due dates come from, and how tight are they?
• Do you need a provable optimum, or a high-quality solution within a time budget — and what is that budget per instance?
• Is a Gurobi license available, or must everything run on open-source tools?
• Will results be benchmarked against Taillard/VRF best-known values for a paper?

Related Skills

• job-shop-scheduling — when jobs have machine-specific routes, or you want the CP-SAT interval model for shop problems with side constraints.
• iterated-local-search — for the ILS framework behind IG: perturbation strength, acceptance criteria, and restart strategies.
• biased-random-key-genetic-algorithm — population-based alternative for the PFSP using a sort decoder over random keys.
• multi-objective-optimization — Pareto methods, epsilon-constraint, NSGA-II, and hypervolume for the makespan-plus-tardiness extension.
• warm-starts-and-initial-solutions — feeding NEH into MIP starts and seeding metaheuristic populations with construction heuristics.