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

```python
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  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 \sumk x{jk} = 1,\ \sumj x{jk} = 1;\qquad c{ik} \ge c{i-1,k} + \textstyle\sumj p{ij} x{jk};\qquad c{ik} \ge c{i,k-1} + \textstyle\sumj 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_co

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.