# Matheuristics > When the user wants to hybridize a MIP solver with heuristic search — fix-and-optimize, relax-and-fix, MIP-based destroy-and-repair (LNS with exact repair), local branching, or solution polishing — including budgeting solver calls inside the loop. Also use when the user mentions "matheuristic," "fix-and-optimize," "relax-and-fix," "local branching," "MIP heuristic hybrid," "proximity search," or when the full MIP stalls while sub-MIPs with most binaries fixed solve in seconds. For heuristic destroy-repair, see large-neighborhood-search; for MIP starts and variable hints, see warm-starts-and-initial-solutions. - Skill: `hajibabaie/matheuristics` (Agent Skill) - Install (CLI): `npx skillmds@latest add hajibabaie/matheuristics` - Raw SKILL.md: https://api.skillmd.com/api/skills/hajibabaie/matheuristics/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Coding & Dev Tools - Author: hajibabaie (https://skillmd.com/u/hajibabaie) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/hajibabaie/matheuristics --- # Matheuristics You are an expert in matheuristics — model-based heuristics that embed an exact MIP solver inside a heuristic search loop. This skill covers fix-and-optimize, relax-and-fix, MIP-based destroy-and-repair (LNS with exact repair), local branching, proximity search and polishing, and the discipline of budgeting solver calls inside the loop. The reference treatment is Maniezzo, Boschetti & Stützle (2021), "Matheuristics: Algorithms and Implementations"; for routing-flavored variants see Archetti & Speranza (2014), "A survey on matheuristics for routing problems." Use the framework below to take a user from "the full MIP stalls at a 5% gap" to a calibrated hybrid that beats both the plain solver and a plain heuristic at equal wall clock. ## Initial Assessment Establish these facts before writing any hybrid code: - **Full-MIP baseline.** Run the complete model with the whole time budget first. Record incumbent, bound, and gap over time. If the solver reaches an acceptable gap, stop — a matheuristic only earns its complexity when the full model stalls. The baseline is also the honesty check every result must be compared against at equal wall clock. - **Where the difficulty lives.** Does the solver struggle to find good incumbents (weak primal side) or to move the bound (weak dual side)? Matheuristics attack the primal side only; if the bound is the problem, look at formulation tightening and cuts instead. - **Decision core.** Which variables are the combinatorial "deciders" (usually binaries: setups, assignments, openings) and which are followers (continuous quantities that an LP determines once the binaries are set)? Fixing schemes fix only the deciders; followers always re-optimize. - **Decomposition dimension.** Is there a natural axis to slice the binaries — time periods, items, machines, regions, vehicles? Fix-and-optimize needs one; local branching and proximity search do not. - **Feasible start.** Does a feasible solution exist already (from a heuristic, from the solver's truncated run, from last week's plan)? If not, construction must be part of the method — relax-and-fix or a truncated MIP run. - **Fixing feasibility.** If part of an incumbent is fixed, does the subproblem stay feasible? Hard capacities, time windows, and inventory balances can make fixings dead-end. If so, plan soft feasibility: penalized slack/overtime variables in the model, kept expensive enough to be driven to zero. - **Sub-MIP sizing.** How many free binaries solve to optimality in 1–5 seconds on this model? Measure it: solve a few random windows of increasing size and log status and runtime. This number drives every window/radius parameter. - **Total budget and split.** How much wall clock per instance, and how should it split between construction, improvement, and a final polishing phase? Count solver calls: budget B with per-call limit tau gives roughly B/tau calls; window schemes must fit. - **Persistent model.** Can one model object live in memory for the whole run so fixing happens through variable bounds? Rebuilding the model per iteration is the most common self-inflicted slowdown. - **Objective structure.** Single objective or lexicographic? Penalty terms already present? Acceptance tests and cutoffs need a single comparable scalar. - **Solver features available.** Gurobi-class solvers expose Cutoff, MIPFocus, Start values, solution pools, and callbacks — all essential here. With CBC/HiGHS the same patterns work but per-call budgets must grow; see the library table. - **Comparison protocol.** Instance set, seeds per instance, and equal-budget reporting against the full MIP — fix these before tuning, exactly as for any metaheuristic. ## Matheuristic Anatomy ### The defining loop All matheuristics in this skill instantiate one template. Given the MIP $\min\{c^\top x + d^\top z : Ax + Bz \ge b,\ x \in \{0,1\}^n,\ z \ge 0\}$ with binaries $x$ and followers $z$, the outer heuristic repeatedly chooses a restriction $N(\bar{x})$ around the incumbent $\bar{x}$ and asks the solver for $$ x^{r+1} \in \arg\min \{ c^\top x + d^\top z : (x,z) \text{ feasible},\ x \in N(x^r) \}, $$ with $N(x^r)$ small enough that the sub-MIP solves within a per-call limit $\tau$, and large enough that it contains solutions the incumbent cannot reach by small heuristic moves. The five standard choices of $N$: | Family | Outer loop chooses | Sub-MIP solves | First reach for | |---|---|---|---| | Fix-and-optimize | a window $W$ of binaries to free | full model with $x_j = \bar{x}_j$ for $j \notin W$ | rolling structures: lot sizing, rostering, multi-period plans | | Relax-and-fix | a stage order (construction, no incumbent needed) | integer stage, LP-relaxed future, fixed past | building the first feasible solution | | MIP-based destroy-repair | a randomized destroy set plus acceptance | exact reinsertion of the destroyed part | tight constraints where greedy repair dead-ends | | Local branching | a Hamming-ball radius $k$ and recentering | model intersected with the ball around $\bar{x}$ | generic binary MIPs with no visible decomposition | | Proximity search / polishing | an improvement threshold $\theta$ | min distance to $\bar{x}$ subject to improving by $\theta$ | refining a good incumbent late in the run | Close relatives worth knowing by name: kernel search (Angelelli, Mansini & Speranza 2010) builds the window from LP information — solve the LP, take variables with nonzero value or small reduced cost as the kernel, then scan "buckets" of the rest; the corridor method (Sniedovich & Voß 2006) imposes a corridor constraint around the incumbent for an exact method other than B&B; POPMUSIC (Taillard & Voß 2002) is fix-and-optimize with parts defined by solution structure. Inside the solver, RINS (Danna, Rothberg & Le Pape 2005, "Exploring relaxation induced neighborhoods to improve MIP solutions") fixes the variables on which the node LP and the incumbent agree and solves the rest — the same idea used automatically at B&B nodes — and solution polishing (Rothberg 2007, "An evolutionary algorithm for polishing mixed integer programming solutions") runs combination/mutation sub-MIPs over a pool of incumbents. ### The two fixing schemes **Fix-and-optimize** (also "exchange"; Pochet & Wolsey 2006, "Production Planning by Mixed Integer Programming"; Helber & Sahling 2010 for the multi-level CLSP) improves a feasible $\bar{x}$. For a window $W$: $$ N_{\text{FO}}(\bar{x}, W) = \{ x \in \{0,1\}^n : x_j = \bar{x}_j \ \ \forall j \notin W \}. $$ Each sub-MIP searches $2^{|W|}$ assignments exactly while the followers $z$ re-optimize globally, so the objective is always the true objective. Windows sweep a decomposition dimension with overlap; passes repeat until a full pass yields no improvement. Monotone descent, no acceptance test needed. **Relax-and-fix** (Pochet & Wolsey 2006; Stadtler 2003 for rolling lot-sizing windows) constructs a solution when none exists. Order the binaries into stages $S_1, \dots, S_m$ (e.g., by time). Stage $s$ solves a model where $S_s$ is integer, $S_1, \dots, S_{s-1}$ are fixed to earlier results, and $S_{s+1}, \dots, S_m$ are LP-relaxed to $[0,1]$. The relaxed future gives lookahead that pure greedy construction lacks. The stage-1 objective is a valid lower bound on the full MIP (it is a relaxation); later stages are not bounds, because fixings have entered. Failure mode: a stage can be infeasible given earlier fixings — counter with stage overlap (re-open the previous stage) or penalized slack in the model. ### The two distance schemes **Local branching** (Fischetti & Lodi 2003, "Local branching") needs no decomposition. Around incumbent $\bar{x}$ with support $S = \{j : \bar{x}_j = 1\}$, the Hamming ball of radius $k$ is the single linear constraint $$ \Delta(x, \bar{x}) = \sum_{j \in S} (1 - x_j) + \sum_{j \notin S} x_j \ \le\ k . $$ Solve the model plus the ball with a node time limit and a cutoff just below $f(\bar{x})$, so only strict improvements are returned. Manage $k$: improvement found — recenter and reset $k$; time-out without improvement — halve $k$ (the ball was too hard to search); ball proven empty — add the reverse cut $\Delta(x,\bar{x}) \ge k+1$ (a valid "tabu" region, since the objective threshold only ever decreases) and enlarge $k$. When the cardinality of the support is constant (assignment/partitioning structure, $\sum_i x_{ij} = 1$), the asymmetric form $\sum_{j \in S}(1-x_j) \le k$ is equivalent up to a factor 2 and halves the constraint length. **Proximity search** (Fischetti & Monaci 2014, "Proximity search for 0-1 mixed-integer convex programming") flips the roles of objective and constraint: add the improvement cut $c^\top x \le c^\top \bar{x} - \theta$, replace the objective by $\min \Delta(x, \bar{x})$, and stop at the first feasible point. The solver's heuristics are very good at "find any solution near here," so each step is cheap; iterate with recentering until $\theta$ cannot be met. ### Decision guidance - **Full MIP first.** Always. It is the baseline, and its truncated run is a free initial incumbent. - **Natural axis + feasible start → fix-and-optimize.** The workhorse for lot sizing, rostering, and multi-period planning. - **No start → relax-and-fix** to construct, then fix-and-optimize to improve. The classic pairing. - **No visible decomposition → local branching** (early, when improvements are easy) and **proximity search** (late, to squeeze a good incumbent). Both treat the model as a black box. - **Destroy logic known, repair hard → MIP-based destroy-repair.** When greedy repair keeps producing infeasible or poor completions, let the solver repair exactly; see large-neighborhood-search for the destroy side. - **Sub-MIPs barely faster than the full model → stop.** If freeing 10% of the binaries already takes half the full-model solve time, the model has little decomposable slack; invest in the formulation instead. ### Parameter guidance | Parameter | Typical range | Increasing it buys | At the cost of | |---|---|---|---| | Window size $|W|$ | 5–25% of binaries; sized so ~80% of sub-MIPs prove optimality within $\tau$ | larger exact moves per solve | longer sub-solves, fewer windows per budget | | Window overlap | 25–50% of window width | improvements that straddle window borders | more solves per pass | | Per-call limit $\tau$ | 1–30 s | more windows searched to proof | fewer calls in the budget | | Sub-MIP gap limit | 0.1–1% (MIPGap per call) | no tail-chasing inside windows | slightly worse window optima | | Passes $P$ | 2–6, stop on a no-improvement pass | extra polishing | diminishing returns | | Relax-and-fix stage width | 2–6 periods (or one resource block) | better lookahead, fewer myopic fixings | slower stages | | Local-branching $k_0$ | 5–20 | richer balls, bigger jumps | node time-outs, $k$ collapses | | Node time limit (LB) | 5–60 s | balls searched to proof, reverse cuts valid | slower recentering | | Proximity $\theta$ | 0.1–1% of incumbent objective | meaningful step per solve | steps prove infeasible sooner | ### Budgeting solver calls The budget arithmetic is unforgiving: passes × windows × $\tau$ must fit the wall clock, with construction and the equal-budget baseline on top. Three rules keep it honest. First, log the status of every sub-solve; if most hit TIME_LIMIT, shrink the window or raise $\tau$ — exact repair that never proves anything is just an expensive heuristic. Second, set `MIPFocus=1` and a per-call `MIPGap` of 0.5–1%: sub-MIPs exist to find good solutions fast, not to prove tight bounds. Third, warm-start every call from the incumbent (`Start` values on the freed variables); the solver then begins with a feasible point and every call is a pure improvement attempt. ## Reusable Engines The two engines below are problem-independent: they see the model through a `dict` mapping a hashable key to its binary variable, and they fix exclusively through variable bounds on one persistent model — the model is never rebuilt. ```text FIX-AND-OPTIMIZE(model M, binaries B, feasible start x̄, windows, passes P, τ) fix every binary to x̄ (by bounds); solve M // continuous completion: f* <- objective // true cost of the start for pass p = 1..P: for each window W in windows(p): // W ⊆ B is the free set unfix W; MIP-start W from x̄ solve M with time limit τ // exact inside W if objective < f* - ε: f* <- objective; x̄|W <- solution re-fix W to x̄ if the pass improved nothing: stop unfix everything; return x̄, f* ``` ```python import math from dataclasses import dataclass, field from typing import Callable, Hashable, Iterable import gurobipy as gp from gurobipy import GRB Key = Hashable @dataclass class FixOptResult: """Incumbent, objective, and per-solve trace of a fix-and-optimize run.""" incumbent: dict[Key, int] objective: float history: list[tuple[int, float, float]] = field(default_factory=list) def solve_sub_mip(model: gp.Model, time_limit: float) -> float: """Solve the current sub-MIP; return its best objective, or +inf if none.""" model.Params.TimeLimit = time_limit model.optimize() ok = model.Status in (GRB.OPTIMAL, GRB.TIME_LIMIT) and model.SolCount > 0 return model.ObjVal if ok else math.inf def fix_and_optimize( model: gp.Model, binaries: dict[Key, gp.Var], start: dict[Key, int], windows: Callable[[int], Iterable[set[Key]]], passes: int = 4, time_per_solve: float = 5.0, tol: float = 1e-6, ) -> FixOptResult: """Improve `start` by re-optimizing one window of binaries at a time. `model` is the full minimization model. `binaries` maps a key to its binary variable; continuous variables are never fixed and re-optimize in every sub-MIP, so objectives are always true objectives. `windows(pass_idx)` yields the key sets left free in each sub-MIP of that pass. Fixing happens through bounds on one persistent model. """ model.Params.OutputFlag = 0 model.Params.MIPFocus = 1 # sub-MIPs exist to find solutions fast incumbent = {k: int(round(v)) for k, v in start.items()} for k, var in binaries.items(): # evaluate the start: everything fixed var.LB = var.UB = incumbent[k] best = solve_sub_mip(model, time_per_solve) if math.isinf(best): raise ValueError(f"start solution infeasible (status {model.Status})") history = [(0, best, model.Runtime)] n_solve = 0 for p in range(passes): improved_in_pass = False for window in windows(p): n_solve += 1 for k in window: # free the window, warm-start it binaries[k].LB, binaries[k].UB = 0.0, 1.0 binaries[k].Start = incumbent[k] obj = solve_sub_mip(model, time_per_solve) if obj < best - tol: best = obj incumbent.update({k: int(round(binaries[k].X)) for k in window}) improved_in_pass = True for k in window: # re-fix to the (possibly new) incumbent binaries[k].LB = binaries[k].UB = incumbent[k] history.append((n_solve, best, model.Runtime)) if not improved_in_pass: # a full pass changed nothing: stop break for var in binaries.values(): # hand back the unfixed model var.LB, var.UB = 0.0, 1.0 return FixOptResult(incumbent, best, history) ``` The companion constructor builds the start that fix-and-optimize needs. Stages are processed in order; each stage is integer, the past is fixed, the future stays LP-relaxed: ```python from typing import Hashable, Sequence import gurobipy as gp from gurobipy import GRB Key = Hashable def relax_and_fix( model: gp.Model, binaries: dict[Key, gp.Var], stages: Sequence[set[Key]], time_per_solve: float = 10.0, overlap: int = 0, ) -> dict[Key, int]: """Construct a feasible solution with a rolling integer window. Stage s solves a MIP in which stages[s] (plus the previous `overlap` stages, re-opened) is binary, later stages are LP-relaxed to [0, 1], and earlier stages are fixed to their solved values. Returns the constructed binary solution; raises RuntimeError if a stage dead-ends (then: increase `overlap`, widen stages, or add penalized slack to the model so fixings cannot cut off feasibility). """ model.Params.OutputFlag = 0 model.Params.TimeLimit = time_per_solve for var in binaries.values(): # everything starts relaxed var.VType = GRB.CONTINUOUS var.LB, var.UB = 0.0, 1.0 fixed: dict[Key, int] = {} for s, stage in enumerate(stages): integer_now = set(stage) for back in range(1, overlap + 1): # re-open recent stages if s - back >= 0: integer_now |= stages[s - back] for k in integer_now: binaries[k].VType = GRB.BINARY if k in fixed: # unfix re-opened keys binaries[k].LB, binaries[k].UB = 0.0, 1.0 model.optimize() if model.SolCount == 0: raise RuntimeError( f"relax-and-fix dead-ended at stage {s} (status {model.Status})") for k in integer_now: fixed[k] = int(round(binaries[k].X)) binaries[k].LB = binaries[k].UB = fixed[k] for var in binaries.values(): # restore the full MIP var.VType = GRB.BINARY var.LB, var.UB = 0.0, 1.0 return fixed ``` ## Worked Application 1: Lot Sizing by Relax-and-Fix + Fix-and-Optimize The capacitated lot-sizing problem (CLSP) is the canonical fix-and-optimize testbed (Trigeiro, Thomas & McClain 1989, "Capacitated lot sizing with setup times"). Items $i$, periods $t$, demand $d_{it}$, setup cost $f_i$ and setup time $st_i$, holding cost $h_i$, capacity $C_t$, overtime $o_t$ at unit cost $c^o$: $$ \min \sum_{i,t} \big( f_i\, y_{it} + h_i\, s_{it} \big) + c^o \sum_t o_t $$ subject to inventory balance $s_{i,t-1} + x_{it} = d_{it} + s_{it}$, the setup link $x_{it} \le M_{it} y_{it}$ with the tight big-M $M_{it} = \sum_{\tau \ge t} d_{i\tau}$, and capacity $\sum_i (x_{it} + st_i\, y_{it}) \le C_t + o_t$. The penalized overtime keeps every sub-MIP feasible no matter what the outer loop fixes — the standard soft-feasibility device in matheuristics. The deciders are the setups $y$; quantities $x, s, o$ are followers. For valid inequalities and reformulations that strengthen this model, see lot-sizing. ```python from dataclasses import dataclass import gurobipy as gp import numpy as np from gurobipy import GRB @dataclass(frozen=True) class ClspInstance: """Capacitated lot sizing: items i = 0..n-1, periods t = 0..T-1.""" demand: np.ndarray # (n, T) setup_cost: np.ndarray # (n,) holding_cost: np.ndarray # (n,) setup_time: np.ndarray # (n,) capacity: np.ndarray # (T,) overtime_cost: float def random_clsp(n_items: int, n_periods: int, seed: int, utilization: float = 0.85) -> ClspInstance: """Trigeiro-style instance; tight capacity, feasibility via overtime.""" rng = np.random.default_rng(seed) demand = rng.integers(40, 141, size=(n_items, n_periods)).astype(float) setup_cost = rng.uniform(300.0, 900.0, n_items) holding_cost = rng.uniform(1.0, 4.0, n_items) setup_time = rng.uniform(10.0, 50.0, n_items) load = demand.sum(axis=0) + setup_time.sum() capacity = np.full(n_periods, utilization * load.mean()) return ClspInstance(demand, setup_cost, holding_cost, setup_time, capacity, overtime_cost=25.0) def build_clsp(inst: ClspInstance) -> tuple[gp.Model, dict[tuple[int, int], gp.Var]]: """CLSP MIP. Returns the model and its setup binaries y[i, t].""" n, T = inst.demand.shape m = gp.Model("clsp") x = m.addVars(n, T, name="x") # production s = m.addVars(n, T, name="s") # end inventory o = m.addVars(T, name="o") # overtime (soft slack) y = m.addVars(n, T, vtype=GRB.BINARY, name="y") # setups: the deciders m.setObjective( gp.quicksum(inst.setup_cost[i] * y[i, t] + inst.holding_cost[i] * s[i, t] for i in range(n) for t in range(T)) + inst.overtime_cost * o.sum(), GRB.MINIMIZE) m.addConstrs( ((s[i, t - 1] if t > 0 else 0.0) + x[i, t] == inst.demand[i, t] + s[i, t] for i in range(n) for t in range(T)), name="balance") rem = inst.demand[:, ::-1].cumsum(axis=1)[:, ::-1] # tight big-M m.addConstrs((x[i, t] <= rem[i, t] * y[i, t] for i in range(n) for t in range(T)), name="setup_link") m.addConstrs( (gp.quicksum(x[i, t] + inst.setup_time[i] * y[i, t] for i in range(n)) <= inst.capacity[t] + o[t] for t in range(T)), name="capacity") m.update() return m, {(i, t): y[i, t] for i in range(n) for t in range(T)} ``` The run constructs with relax-and-fix over 6-period time stages (overlap 1), then improves with fix-and-optimize, alternating time windows and item windows between passes — alternating the decomposition dimension is what lets the method escape the blind spots of either slicing: ```python import time def time_windows(n_items: int, n_periods: int, width: int, step: int) -> list[set[tuple[int, int]]]: """Overlapping windows over periods: all items, t in [t0, t0 + width).""" return [{(i, t) for i in range(n_items) for t in range(t0, min(t0 + width, n_periods))} for t0 in range(0, max(n_periods - width, 0) + 1, step)] def item_windows(n_items: int, n_periods: int, block: int) -> list[set[tuple[int, int]]]: """Windows over items: all periods for `block` items at a time.""" return [{(i, t) for i in range(i0, min(i0 + block, n_items)) for t in range(n_periods)} for i0 in range(0, n_items, block)] inst = random_clsp(n_items=30, n_periods=60, seed=7) # 1800 setup binaries model, y = build_clsp(inst) n_items, n_periods = inst.demand.shape t0 = time.time() stages = [{(i, t) for i in range(n_items) for t in range(s0, min(s0 + 6, n_periods))} for s0 in range(0, n_periods, 6)] start = relax_and_fix(model, y, stages, time_per_solve=5.0, overlap=1) def clsp_windows(p: int) -> list[set[tuple[int, int]]]: """Alternate the decomposition dimension between passes.""" if p % 2 == 0: return time_windows(n_items, n_periods, width=6, step=3) return item_windows(n_items, n_periods, block=5) result = fix_and_optimize(model, y, start, clsp_windows, passes=4, time_per_solve=3.0) wall = time.time() - t0 full, _ = build_clsp(inst) # honest equal-budget comparison full.Params.OutputFlag = 0 full.Params.TimeLimit = wall full.optimize() assert full.SolCount > 0 print(f"relax-and-fix start {result.history[0][1]:.0f} -> " f"fix-and-optimize {result.objective:.0f} " f"({len(result.history) - 1} sub-MIPs, {wall:.0f}s)") print(f"full MIP at {wall:.0f}s: obj {full.ObjVal:.0f}, " f"bound {full.ObjBound:.0f}") # Expected (Gurobi 13; time-limit-driven, so machine-dependent): relax-and- # fix start ~814000; fix-and-optimize ~785000 after ~25 sub-MIPs, ~125s in # total. The full MIP given the same wall clock returns ~787000 over bound # ~776000 -- the hybrid wins at equal budget, with margin growing in size. ``` Two design points carry the example. First, the per-solve trace in `result.history` shows which windows produced the improvements; in CLSP almost all gains come from time windows around capacity-tight periods, which is the signal to add more overlap exactly there. Second, the equal-budget full-MIP line is not decoration — without it, a matheuristic result is unreviewable. ## Worked Application 2: Generalized Assignment by Local Branching The generalized assignment problem (GAP) assigns each job $j$ to exactly one agent $i$, minimizing $\sum_{ij} c_{ij} x_{ij}$ subject to $\sum_i x_{ij} = 1$ and agent budgets $\sum_j a_{ij} x_{ij} \le b_i$. There is no natural window dimension — capacity couples every job on an agent — so local branching is the right matheuristic. Because $\sum_i x_{ij} = 1$ fixes the support size at $n$, the asymmetric ball $\sum_{(i,j): \bar{x}_{ij}=1} (1 - x_{ij}) \le k$ reads "at most $k$ jobs change agent." ```text LOCAL-BRANCHING(model M, binaries B, k0, node time τ, total budget T) x̄ <- incumbent from a truncated run of M // any feasible point works k <- k0 while budget remains and 1 <= k <= |support(x̄)|: add ball Δ(x, x̄) <= k; set Cutoff just below f(x̄) solve M for at most τ seconds improved -> recenter at the new incumbent; k <- k0 time-out, none -> k <- k/2 // ball too large to search proven empty -> add reverse cut Δ(x, x̄) >= k+1; k <- k + k0/2 remove the ball return the best incumbent found ``` ```python from dataclasses import dataclass import gurobipy as gp import numpy as np from gurobipy import GRB @dataclass(frozen=True) class GapInstance: """Generalized assignment: agents i = 0..m-1, jobs j = 0..n-1.""" cost: np.ndarray # (m, n) resource: np.ndarray # (m, n) budget: np.ndarray # (m,) def random_gap(n_agents: int, n_jobs: int, seed: int, tightness: float = 0.8) -> GapInstance: """OR-Library class-D style: cost anti-correlated with resource (hard).""" rng = np.random.default_rng(seed) resource = rng.integers(1, 101, size=(n_agents, n_jobs)).astype(float) cost = 111.0 - resource + rng.integers(-10, 11, size=(n_agents, n_jobs)) budget = tightness * resource.sum(axis=1) / n_agents return GapInstance(cost, resource, budget) def build_gap(inst: GapInstance) -> tuple[gp.Model, dict[tuple[int, int], gp.Var]]: """GAP MIP: every job to exactly one agent, within agent budgets.""" m_agents, n_jobs = inst.cost.shape m = gp.Model("gap") x = m.addVars(m_agents, n_jobs, vtype=GRB.BINARY, name="x") m.setObjective( gp.quicksum(inst.cost[i, j] * x[i, j] for i in range(m_agents) for j in range(n_jobs)), GRB.MINIMIZE) m.addConstrs((x.sum("*", j) == 1 for j in range(n_jobs)), name="assign") m.addConstrs( (gp.quicksum(inst.resource[i, j] * x[i, j] for j in range(n_jobs)) <= inst.budget[i] for i in range(m_agents)), name="budget") m.update() return m, {(i, j): x[i, j] for i in range(m_agents) for j in range(n_jobs)} ``` The loop uses the solver's `Cutoff` parameter so a sub-solve can only return strict improvements; a proven-empty ball is then status `CUTOFF` (or `INFEASIBLE` once reverse cuts accumulate), which is a proof, not an error: ```python from typing import Hashable import gurobipy as gp from gurobipy import GRB Key = Hashable def local_branching( model: gp.Model, binaries: dict[Key, gp.Var], k0: int = 15, node_time: float = 5.0, total_time: float = 60.0, ) -> tuple[dict[Key, int], float, list[tuple[float, float]]]: """Local branching (Fischetti & Lodi 2003) for a minimization MIP. Uses the asymmetric Hamming ball over the incumbent support, valid whenever the number of ones is constant across feasible solutions (assignment / partitioning structure). Returns (best solution, objective, trace of (clock, incumbent objective)).""" model.Params.OutputFlag = 0 model.Params.MIPFocus = 1 model.Params.TimeLimit = min(node_time, total_time) model.optimize() # truncated run -> first incumbent if model.SolCount == 0: raise RuntimeError(f"no initial solution (status {model.Status})") center = {key: int(round(var.X)) for key, var in binaries.items()} best, best_obj = dict(center), model.ObjVal clock = model.Runtime trace = [(clock, best_obj)] k, n_tabu = k0, 0 n_ones = sum(center.values()) while clock < total_time and 1 <= k <= n_ones: support = [binaries[key] for key, v in center.items() if v == 1] ball = model.addConstr(gp.quicksum(1 - v for v in support) <= k, name="lb_ball") model.Params.TimeLimit = min(node_time, total_time - clock) model.Params.Cutoff = best_obj - 1e-6 # improvements only model.optimize() clock += model.Runtime if model.SolCount > 0: # strictly better by construction center = {key: int(round(var.X)) for key, var in binaries.items()} best, best_obj = dict(center), model.ObjVal model.remove(ball) k = k0 # recenter, reset the radius elif model.Status == GRB.TIME_LIMIT: model.remove(ball) k //= 2 # ball too hard to search: shrink else: # CUTOFF / INFEASIBLE: proven empty model.remove(ball) n_tabu += 1 # forbid the explored ball for good model.addConstr(gp.quicksum(1 - v for v in support) >= k + 1, name=f"lb_tabu_{n_tabu}") k += max(1, k0 // 2) # and look further out trace.append((clock, best_obj)) model.Params.Cutoff = GRB.INFINITY # leave the model clean return best, best_obj, trace ``` The reverse cut deserves one line of justification: when a ball around an old center is proven to contain nothing better than the then-incumbent, it can never contain anything better than any later (smaller) incumbent, so the cut stays valid for the rest of the run. Run it on a hard class-D instance, validate independently, and compare at equal budget: ```python import numpy as np gap_inst = random_gap(n_agents=10, n_jobs=100, seed=3) gap_model, x = build_gap(gap_inst) best, obj, trace = local_branching(gap_model, x, k0=15, node_time=5.0, total_time=45.0) # Independent validation -- never trust the loop's own bookkeeping. assign = {j: i for (i, j), v in best.items() if v == 1} loads = np.zeros(gap_inst.cost.shape[0]) for j, i in assign.items(): loads[i] += gap_inst.resource[i, j] assert len(assign) == gap_inst.cost.shape[1] assert np.all(loads <= gap_inst.budget + 1e-6) assert abs(sum(gap_inst.cost[i, j] for j, i in assign.items()) - obj) < 1e-6 ref, _ = build_gap(gap_inst) ref.Params.OutputFlag = 0 ref.Params.TimeLimit = 45.0 ref.optimize() print(f"local branching {trace[0][1]:.0f} -> {obj:.0f}; " f"full MIP 45s: {ref.ObjVal:.0f} (bound {ref.ObjBound:.0f})") # Expected (Gurobi 13; time-limit-driven, so run-dependent): initial # incumbent ~6335-6340; local branching improves it to ~6324-6331, on par # with or better than the full MIP at equal wall clock (bound ~6318). The # trace shows the k-management cycle: improve (k resets), time-out # (k halves), proven-empty ball (reverse cut added, k grows). ``` ## Advanced Techniques ### MIP-based destroy-and-repair Fix-and-optimize with a *randomized* window is exactly LNS with exact repair: destroy = unfix a coherent group of binaries, repair = sub-MIP. The engine below takes the groups (one group per job, customer, period, or machine — the removable element of large-neighborhood-search) and hill-climbs over exact large moves; add Metropolis acceptance or adaptive group-family weights exactly as in large-neighborhood-search when it stalls. On the GAP instance above, freeing 40 of 100 job columns for 2-second repairs improves a truncated-run incumbent by several units within seconds. ```python from typing import Hashable, Sequence import gurobipy as gp import numpy as np from gurobipy import GRB Key = Hashable def mip_lns( model: gp.Model, binaries: dict[Key, gp.Var], start: dict[Key, int], groups: Sequence[Sequence[Key]], q: int, iterations: int = 50, time_per_solve: float = 2.0, seed: int = 0, ) -> tuple[dict[Key, int], float]: """LNS with exact repair: free q random groups of binaries, re-solve. Improvement-only acceptance; the destroyed part is repaired optimally (within the time limit) while everything else stays fixed by bounds.""" rng = np.random.default_rng(seed) model.Params.OutputFlag = 0 model.Params.MIPFocus = 1 model.Params.TimeLimit = time_per_solve incumbent = {k: int(round(v)) for k, v in start.items()} for k, var in binaries.items(): var.LB = var.UB = incumbent[k] model.optimize() if model.SolCount == 0: raise ValueError("start solution infeasible under full fixing") best = model.ObjVal for _ in range(iterations): chosen = rng.choice(len(groups), size=q, replace=False) window = [k for g in chosen for k in groups[g]] for k in window: binaries[k].LB, binaries[k].UB = 0.0, 1.0 binaries[k].Start = incumbent[k] model.optimize() if model.SolCount > 0 and model.ObjVal < best - 1e-6: best = model.ObjVal for k in window: incumbent[k] = int(round(binaries[k].X)) for k in window: binaries[k].LB = binaries[k].UB = incumbent[k] for var in binaries.values(): var.LB, var.UB = 0.0, 1.0 return incumbent, best ``` The destroy set should be *related* — groups that share a constraint (same agent, adjacent periods, nearby customers) — because the sub-MIP can only trade among what is simultaneously free. Random groups are the baseline; upgrade the selection with the relatedness measures from large-neighborhood-search. ### Proximity search and polishing When fix-and-optimize and local branching both stall, proximity search often still moves: asking the solver for *any* solution at least $\theta$ better, as close as possible to the incumbent, exploits exactly what MIP heuristics are best at. `SolutionLimit=1` makes each step cheap. ```python from typing import Hashable import gurobipy as gp from gurobipy import GRB Key = Hashable def proximity_search_step( model: gp.Model, binaries: dict[Key, gp.Var], center: dict[Key, int], center_obj: float, theta: float, time_limit: float = 10.0, ) -> dict[Key, int] | None: """One proximity-search step (Fischetti & Monaci 2014). Adds the improvement cut obj <= center_obj - theta, minimizes the Hamming distance to `center`, stops at the first feasible point. Returns the new center, or None if no improving point was found.""" original = model.getObjective() cutoff = model.addConstr(original <= center_obj - theta, name="ps_cutoff") dist = gp.quicksum((1 - v) if center[k] == 1 else v for k, v in binaries.items()) model.setObjective(dist, GRB.MINIMIZE) model.Params.OutputFlag = 0 model.Params.TimeLimit = time_limit model.Params.SolutionLimit = 1 # the first improving point suffices model.optimize() new_center = None if model.SolCount > 0: new_center = {k: int(round(v.X)) for k, v in binaries.items()} model.remove(cutoff) # restore the original model model.setObjective(original, GRB.MINIMIZE) model.Params.SolutionLimit = GRB.MAXINT return new_center ``` Iterate with recentering until a step returns `None` at the smallest useful $\theta$. This is also the cheap "polishing phase" to schedule for the last 10–20% of any matheuristic budget. Solvers ship their own polishing: Gurobi's `ImproveStartTime`/`ImproveStartGap` switch the B&B to pure solution improvement (Rothberg 2007), which is a zero-code alternative worth benchmarking against. ### Soft fixing and LP-guided windows Hard fixing ($x_j = \bar{x}_j$) can be replaced by *soft fixing*: a local-branching constraint over the would-be-fixed set, allowing the sub-MIP to overrule a few fixings where it pays. Fischetti & Lodi's $\Delta \le k$ over the complement of the window is the standard form. LP information picks better windows than blind sweeps: solve the LP relaxation once per pass and (a) free the variables whose LP value disagrees with the incumbent (the RINS complement — Danna, Rothberg & Le Pape 2005), or (b) free variables with reduced cost near zero, the kernel-search criterion (Angelelli, Mansini & Speranza 2010). Both target exactly the variables the relaxation says are decided wrongly or barely decided. ### Adaptive window selection Window families (time, item, machine, LP-guided, random-related) are operators; their per-family improvement-per-second statistics are rewards. Reuse the ALNS roulette and segment-update machinery from large-neighborhood-search verbatim, with one practical change: normalize rewards by sub-solve runtime, because a 30-second window that saves 100 cost units is worse than five 3-second windows saving 30 each. ### Warm starts and budget shaping across calls Every sub-solve should start from the incumbent: set `Start` on freed variables (fixed variables are determined by bounds). For sub-MIPs that mostly prove optimality early, cap effort with a per-call `MIPGap` (0.5–1%) in addition to `TimeLimit`, and consider `BestObjStop` to return as soon as any improvement exists. Across the run, shape budgets like a cooling schedule: short calls and small windows early (many cheap improvements), longer calls and larger windows late (rare deep improvements), proximity-search polishing last. The mechanics of starts, hints, and partial starts are covered in warm-starts-and-initial-solutions; callback-level control (aborting a sub-solve once improvement stalls) in gurobi-advanced-features. ## Practical Challenges **Sub-MIPs hit the time limit on every window.** The window is too large for the model's difficulty, or the sub-MIP inherits the full model's weak relaxation. Shrink the window until ~80% of calls prove optimality; set a per-call `MIPGap` so near-proofs terminate; check that fixing is done by bounds (presolve then removes the fixed variables — fixing by added equality constraints is much weaker). If even tiny windows are slow, the formulation is the problem, not the matheuristic. **Relax-and-fix dead-ends or returns absurdly expensive plans.** A stage fixed itself into a corner. Add penalized slack variables (overtime, demand shortfall at high cost) so no fixing can cut off feasibility; increase stage overlap so the previous stage can be revised; widen stages so the integer window sees further. If the construction is feasible but poor, the LP-relaxed future is too optimistic — shorten stages or keep one extra stage integer. **Fix-and-optimize stalls after the first pass.** Single-dimension windows converge to a "window-optimal" fixed point quickly. Alternate decomposition dimensions between passes, add 30–50% overlap, randomize window order, and grow window width by 50% whenever a pass fails — escalation is cheap because the pass would otherwise terminate the run anyway. Soft fixing (see Advanced Techniques) is the next escalation step. **The first local-branching ball already times out.** $k_0$ is too large or the node budget too small; the ball sub-MIP is nearly as hard as the full model. Halve $k_0$, raise the node limit, use the asymmetric ball when support size is constant, and always warm-start the ball with the center. If the model is hard even at $k = 5$, switch to proximity search, which only needs the first feasible point. **Gurobi reports CUTOFF and the loop treats it as a failure.** With `Cutoff` set, "no solution" is the *expected* outcome of a proven-empty ball, status `GRB.CUTOFF`. Handle the triple (improved / `TIME_LIMIT` without solution / prov …(truncated)