diff --git a/CHANGELOG.md b/CHANGELOG.md index ae37c3270..c56c027c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,8 @@ ## Unreleased ### Added +- Added example `heur_lpt.py`: warm-starting a scheduling MIP with a `Heur` plugin that builds its solution in the original space +- Added example `scheduling_logical.py`: the same scheduling problem modelled with and/or/indicator, disjunction and cardinality constraints - Added methods: `getNNodesLeft()`, `getNRuns()`, `getNReoptRuns()`, `addNNodes()` with tests - Added `addConsCumulative()` for SCIP cumulative constraints (#1222) - `Expr` and `GenExpr` support `__pos__` magic method like `+Expr` or `+GenExpr` diff --git a/examples/finished/heur_lpt.py b/examples/finished/heur_lpt.py new file mode 100644 index 000000000..30e66c01b --- /dev/null +++ b/examples/finished/heur_lpt.py @@ -0,0 +1,157 @@ +""" +Parallel machine scheduling with release dates (Pm|r_j|C_max), solved with a +big-M formulation that is warm-started by a custom heuristic. + +The heuristic is a Heur plugin that runs once before the root node. It builds +an LPT (longest processing time first) list schedule and hands it to SCIP as +an incumbent. +""" + +from pyscipopt import Model, Heur, SCIP_RESULT, SCIP_HEURTIMING, quicksum + + +def lpt_schedule(p, r, m): + """ + LPT list scheduling: whenever a machine becomes free, start the longest + released job on it. Returns {job: (machine, start)} and the makespan. + """ + unscheduled = set(range(len(p))) + free_at = [0] * m + schedule = {} + + while unscheduled: + machine = min(range(m), key=lambda k: free_at[k]) + now = free_at[machine] + + released = [j for j in unscheduled if r[j] <= now] + if not released: + now = min(r[j] for j in unscheduled) + released = [j for j in unscheduled if r[j] <= now] + + job = max(released, key=lambda j: p[j]) + schedule[job] = (machine, now) + free_at[machine] = now + p[job] + unscheduled.remove(job) + + return schedule, max(free_at) + + +class LPTHeur(Heur): + """Offers the LPT schedule to SCIP as a solution.""" + + def __init__(self, p, r, m, start, assign, before, makespan): + super().__init__() + self.p = p + self.r = r + self.m = m + self.start = start + self.assign = assign + self.before = before + self.makespan = makespan + self.done = False + + def heurexec(self, heurtiming, nodeinfeasible): + # run once; SCIP processes the root again after a restart + if self.done: + return {"result": SCIP_RESULT.DIDNOTRUN} + self.done = True + + schedule, cmax = lpt_schedule(self.p, self.r, self.m) + + # Build the solution in the original space. Presolving may have fixed + # some variables (symmetry handling does so here), and setting a + # conflicting value on a fixed variable of a transformed solution fails. + sol = self.model.createOrigSol(self) + + sol[self.makespan] = cmax + for j, (machine, st) in schedule.items(): + sol[self.start[j]] = st + for k in range(self.m): + sol[self.assign[j, k]] = 1 if k == machine else 0 + + for i, (machine_i, st_i) in schedule.items(): + for j, (machine_j, st_j) in schedule.items(): + if i != j: + same_machine = machine_i == machine_j + sol[self.before[i, j]] = 1 if same_machine and st_i < st_j else 0 + + accepted = self.model.trySol(sol) + print(f"LPT heuristic: makespan {cmax} {'accepted' if accepted else 'rejected'}") + + if accepted: + return {"result": SCIP_RESULT.FOUNDSOL} + return {"result": SCIP_RESULT.DIDNOTFIND} + + +def build_model(p, r, m): + """ + Disjunctive big-M model. + + start[j] - start time of job j + assign[j, k] - 1 if job j runs on machine k + before[i, j] - 1 if job i finishes before job j starts + makespan - completion time of the last job + """ + jobs = range(len(p)) + machines = range(m) + horizon = max(r) + sum(p) # no job needs to start later than this + + model = Model("Pm|r_j|C_max") + + start = {j: model.addVar(vtype="C", lb=r[j], ub=horizon - p[j], name=f"start_{j}") for j in jobs} + assign = {(j, k): model.addVar(vtype="B", name=f"assign_{j}_{k}") for j in jobs for k in machines} + before = {(i, j): model.addVar(vtype="B", name=f"before_{i}_{j}") for i in jobs for j in jobs if i != j} + makespan = model.addVar(vtype="C", ub=horizon, name="makespan") + + for j in jobs: + model.addCons(quicksum(assign[j, k] for k in machines) == 1) + model.addCons(start[j] + p[j] <= makespan) + + for i in jobs: + for j in jobs: + if i != j: + model.addCons(start[i] + p[i] <= start[j] + horizon * (1 - before[i, j])) + + # jobs on the same machine have to be sequenced + for i in jobs: + for j in jobs: + if i < j: + for k in machines: + model.addCons(assign[i, k] + assign[j, k] - 1 <= before[i, j] + before[j, i]) + + model.setObjective(makespan, "minimize") + + return model, start, assign, before, makespan + + +def print_schedule(p, schedule, m): + for k in range(m): + on_k = sorted((st, j) for j, (machine, st) in schedule.items() if machine == k) + jobs = " ".join(f"job {j} [{st}, {st + p[j]})" for st, j in on_k) + print(f" machine {k}: {jobs}") + + +if __name__ == "__main__": + p = [2, 2, 2, 9, 8, 7, 4, 5, 3, 6] + r = [0, 0, 0, 1, 1, 1, 3, 5, 2, 4] + m = 3 + + schedule, cmax = lpt_schedule(p, r, m) + print(f"LPT schedule, makespan {cmax}:") + print_schedule(p, schedule, m) + + model, start, assign, before, makespan = build_model(p, r, m) + + heur = LPTHeur(p, r, m, start, assign, before, makespan) + model.includeHeur(heur, "lpt", "LPT list scheduling warm start", "L", + freq=0, timingmask=SCIP_HEURTIMING.BEFORENODE) # freq=0: root node only + + model.optimize() + + print(f"\nSCIP status: {model.getStatus()}") + print(f"optimal schedule, makespan {model.getObjVal():g}:") + best = {} + for j in range(len(p)): + machine = next(k for k in range(m) if model.getVal(assign[j, k]) > 0.5) + best[j] = (machine, round(model.getVal(start[j]))) + print_schedule(p, best, m) diff --git a/examples/finished/scheduling_logical.py b/examples/finished/scheduling_logical.py new file mode 100644 index 000000000..29e3a08c4 --- /dev/null +++ b/examples/finished/scheduling_logical.py @@ -0,0 +1,150 @@ +""" +Parallel machine scheduling with release dates (Pm|r_j|C_max), modelled three +ways with SCIP's logical constraints instead of big-M constraints: + + and_or_indicator_model and, or and indicator constraints + disjunction_model disjunction constraints + time_indexed_model cardinality constraints + +The formulations differ in what they contribute to the LP relaxation, which +shows in the number of nodes SCIP needs. +""" + +from pyscipopt import Model, quicksum + + +def base_model(p, r, m, name): + """Everything except the condition that jobs on one machine do not overlap.""" + jobs = range(len(p)) + machines = range(m) + horizon = max(r) + sum(p) # no job needs to start later than this + + model = Model(name) + + start = {j: model.addVar(vtype="C", lb=r[j], ub=horizon - p[j], name=f"start_{j}") for j in jobs} + assign = {(j, k): model.addVar(vtype="B", name=f"assign_{j}_{k}") for j in jobs for k in machines} + makespan = model.addVar(vtype="C", ub=horizon, name="makespan") + + for j in jobs: + # not an xor constraint: xor fixes the parity, so with three machines + # a job could be assigned to all of them + model.addCons(quicksum(assign[j, k] for k in machines) == 1) + model.addCons(start[j] + p[j] <= makespan) + + model.setObjective(makespan, "minimize") + + return model, start, assign, makespan + + +def and_or_indicator_model(p, r, m): + """ + before[i, j] = 1 forces j to start after i is done (indicator). For each + pair and machine: same = assign[i, k] and assign[j, k], ordered = + before[i, j] or before[j, i], and same implies ordered (indicator). + """ + model, start, assign, makespan = base_model(p, r, m, "and-or-indicator") + jobs = range(len(p)) + machines = range(m) + + before = {(i, j): model.addVar(vtype="B", name=f"before_{i}_{j}") for i in jobs for j in jobs if i != j} + + for i in jobs: + for j in jobs: + if i != j: + model.addConsIndicator(start[i] + p[i] <= start[j], binvar=before[i, j], name=f"seq_{i}_{j}") + + for i in jobs: + for j in jobs: + if i < j: + ordered = model.addVar(vtype="B", name=f"ordered_{i}_{j}") + model.addConsOr([before[i, j], before[j, i]], ordered) + for k in machines: + same = model.addVar(vtype="B", name=f"same_{i}_{j}_{k}") + model.addConsAnd([assign[i, k], assign[j, k]], same) + model.addConsIndicator(ordered >= 1, binvar=same, name=f"ordered_if_same_{i}_{j}_{k}") + + return model, start, assign, makespan + + +def disjunction_model(p, r, m): + """ + For each pair and machine: i is not on the machine, or j is not, or i + finishes before j starts, or j finishes before i starts. + + A disjunction is enforced by branching only and adds nothing to the LP + relaxation, so without the load constraint at the end SCIP enumerates + schedules. Try removing it. + """ + model, start, assign, makespan = base_model(p, r, m, "disjunction") + jobs = range(len(p)) + machines = range(m) + + for i in jobs: + for j in jobs: + if i < j: + for k in machines: + model.addConsDisjunction( + [assign[i, k] <= 0, assign[j, k] <= 0, + start[i] + p[i] <= start[j], start[j] + p[j] <= start[i]], + name=f"no_overlap_{i}_{j}_{k}", + ) + + for k in machines: + # total processing time on a machine is a lower bound on the makespan + model.addCons(quicksum(p[j] * assign[j, k] for j in jobs) <= makespan, name=f"load_{k}") + + return model, start, assign, makespan + + +def time_indexed_model(p, r, m): + """ + x[j, k, t] = 1 if job j starts on machine k at time t. A cardinality + constraint per machine and time step allows at most one running job. + """ + model, start, assign, makespan = base_model(p, r, m, "time-indexed") + jobs = range(len(p)) + machines = range(m) + horizon = max(r) + sum(p) + slots = {j: range(r[j], horizon - p[j] + 1) for j in jobs} # possible start times + + x = {} + for j in jobs: + for k in machines: + for t in slots[j]: + x[j, k, t] = model.addVar(vtype="B", name=f"x_{j}_{k}_{t}") + + for j in jobs: + model.addCons(quicksum(x[j, k, t] for k in machines for t in slots[j]) == 1) + model.addCons(start[j] == quicksum(t * x[j, k, t] for k in machines for t in slots[j])) + for k in machines: + model.addCons(assign[j, k] == quicksum(x[j, k, t] for t in slots[j])) + + for k in machines: + for t in range(horizon): + # jobs that would be running on machine k at time t + running = [x[j, k, s] for j in jobs for s in range(t - p[j] + 1, t + 1) if (j, k, s) in x] + if len(running) > 1: + model.addConsCardinality(running, 1, name=f"one_job_{k}_{t}") + + return model, start, assign, makespan + + +def print_schedule(p, m, start, assign, model): + for k in range(m): + on_k = sorted((round(model.getVal(start[j])), j) for j in range(len(p)) if model.getVal(assign[j, k]) > 0.5) + jobs = " ".join(f"job {j} [{st}, {st + p[j]})" for st, j in on_k) + print(f" machine {k}: {jobs}") + + +if __name__ == "__main__": + p = [2, 2, 2, 9, 8, 7, 4, 5, 3, 6] + r = [0, 0, 0, 1, 1, 1, 3, 5, 2, 4] + m = 3 + + for build in (and_or_indicator_model, disjunction_model, time_indexed_model): + model, start, assign, makespan = build(p, r, m) + model.hideOutput() + model.optimize() + print(f"{model.getProbName()}: {model.getStatus()}, makespan {model.getObjVal():g}, " + f"{model.getNNodes()} nodes, {model.getSolvingTime():.2f}s") + print_schedule(p, m, start, assign, model)