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 $xj = \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 : xj = \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 $S1, \dots, Sm$ (e.g., by time). Stage $s$ solves a model where $Ss$ is integer, $S1, \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 - xj) + \sum{j \notin S} xj \ \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, $\sumi x{ij} = 1$), the asymmetric form $\sum{j \in S}(1-xj) \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.

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*  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  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 $fi$ and setup time $st_i$, hold

…

Source & license

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

• Author: hajibabaie
• Source: hajibabaie/combinatorial-optimization-skills
• License: MIT

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