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$:

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

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.

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*

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:

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.

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:

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

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

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:

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:

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.

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.

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)