MILP Modeling with Gurobi

You are an expert in mathematical optimization and the gurobipy API. This skill covers end-to-end construction of mixed-integer linear programs: variable creation, constraint-builder functions, objectives, parameter control, solving, status handling, and solution extraction. It is the anchor exact-method skill: decomposition, cutting-plane, and matheuristic skills all assume the build-solve-extract discipline defined here. Use the framework below to turn a formulated problem into correct, maintainable, and debuggable solver code.

Initial Assessment

Establish these facts before writing any model code:

• Problem class and formulation status. Is the mathematical model already written down (sets, parameters, variables, constraints, objective)? If not, formulate first — see problem-formulation. Never start typing addVar against a vague word problem.
• Instance size. Count variables and constraints as functions of the data (e.g., GAP has m*n binaries, n + m rows). Estimate nonzeros. Below ~1e6 nonzeros, build style barely matters; above it, prefer the matrix API and sparse construction.
• Integrality. Which decisions are truly discrete? Every avoidable integer variable costs branching effort. Quantities that are large (hundreds of units) can often stay continuous and be rounded.
• Hard vs soft constraints. Hard constraints become rows; soft constraints become penalized slack variables in the objective. Confirm the classification with the user before coding — it changes the model skeleton.
• License and solver availability. Gurobi requires a license (free academic licenses exist; the pip install ships a size-limited trial). If there is no license, switch to open-source-solvers before investing in gurobipy-specific code.
• Time budget and quality target. Proven optimality, 1% gap, or "best in 60 seconds"? This sets MIPGap and TimeLimit up front and decides whether a heuristic is the better tool entirely.
• Data format and units. Where do coefficients come from (CSV, database, generated)? Check unit consistency: mixing cents with millions of dollars creates the numerical pathologies described under Practical Challenges.
• Re-solve pattern. One-shot solve, or repeated solves with changing data (rolling horizon, decomposition loop)? Repeated solves favor model modification over rebuild and explicit Seed/Threads control for reproducibility.
• Deliverable. A solution table, a reusable module, or an experiment harness? This determines how much of the Output Format section applies.

Model Anatomy

The object of study is the mixed-integer linear program

$$ \min_{x,,y} ; c^{\top} x + d^{\top} y \quad \text{s.t.} \quad A x + G y \le b, \qquad x \in \mathbb{Z}^{n}{\ge 0}, \quad y \in \mathbb{R}^{m}{\ge 0}. $$

Dropping integrality gives the LP relaxation, solvable in polynomial time; the integrality restriction makes MILP NP-hard, and branch-and-bound explores a worst-case exponential tree. In practice difficulty is governed less by size than by the tightness of the LP relaxation: a 1e6-variable model with a tight relaxation can solve in minutes while a 500-variable model with loose big-M constraints stalls for hours (Wolsey (1998), Integer Programming; Klotz & Newman (2013), Practical guidelines for solving difficult mixed integer linear programs). Formulation strength is covered in integer-programming-techniques; this skill makes sure the model you build is the model you meant.

The gurobipy object model

Object	Created by	Role
gp.Env	implicit, or gp.Env() explicitly	license session; create explicitly in services so it can be disposed cleanly
gp.Model	gp.Model("name")	owns variables, constraints, objective, parameters, attributes
gp.Var	model.addVar(...)	single decision variable with lb, ub, vtype, obj, name
gp.tupledict of Var	model.addVars(...)	indexed variable family; supports .sum(...), .select(...), .items()
gp.LinExpr	arithmetic, gp.quicksum	linear expression; quicksum is the fast accumulator
gp.Constr	model.addConstr(s)(...)	named linear row; query Pi, Slack, RHS after solve
gp.MVar / matrix constraints	model.addMVar, A @ x <= b	numpy-style blocks for large regular structure

Variable types

vtype	Use for	Watch out
GRB.CONTINUOUS	quantities, flows, slacks	default lb is 0.0; free variables need lb=-GRB.INFINITY
GRB.BINARY	on/off, assignment, setup decisions	solution values are floats near 0/1; test with > 0.5, never == 1
GRB.INTEGER	counts, lot sizes	always set a finite, data-derived ub; loose bounds weaken the relaxation
GRB.SEMICONT / GRB.SEMIINT	value is 0 or in [lb, ub]	replaces an explicit binary + big-M link; solver branches on the disjunction

Lifecycle rules

• Lazy updates. gurobipy batches addVar/addConstr calls; the model is synchronized at optimize() or an explicit model.update(). Querying an attribute of a just-added object before an update raises an error — this is the classic first-session surprise.
• Attributes are phase-dependent. var.X, model.ObjVal, constr.Pi exist only after a solve that produced what you are asking for. Gate every read on model.Status and model.SolCount as in the harness below.
• Name everything. Named variables and constraints make model.write("model.lp") readable. Writing the LP file and reading it line by line is the single most effective model-debugging technique: you see the model the solver sees, after your indexing bugs.
• Modify, don't rebuild. constr.RHS, var.Obj, var.LB/UB, and model.chgCoeff support re-solve loops; Gurobi reuses the LP basis or MIP incumbent where possible.

When MILP is the right hammer

• Objective and constraints are linear, or linearizable at acceptable cost (products, absolute values, piecewise terms — see linearization-techniques).
• Proven bounds or optimality certificates matter (publication, money, safety).
• Instance sizes are moderate, or the relaxation is known to be strong.
• Prefer CP-SAT for deeply disjunctive scheduling, and metaheuristics when the evaluation budget is tiny or constraints resist linearization; the triage logic lives in problem-formulation.

A Reusable Build-Solve-Extract Skeleton

Every model in this repository follows one skeleton. Constraint families get one builder function each — this keeps the formulation reviewable (each function maps to one block of math), unit-testable, and reusable across model variants.

BUILD-SOLVE-EXTRACT(data):
1. model <- Model(name); set Params explicitly (TimeLimit, MIPGap, Seed, Threads)
2. vars  <- add_variables(model, data)          # vtype, lb, ub fixed at creation time
3. for each constraint family F in the formulation:
       add_F_constraints(model, vars, data)     # one function per family, named rows
4. set_objective(model, vars, data)
5. model.optimize()
6. case model.Status of
       OPTIMAL                      -> extract ObjVal, var.X, ObjBound, MIPGap
       TIME_LIMIT and SolCount > 0  -> extract incumbent, report remaining gap
       TIME_LIMIT and SolCount == 0 -> report bound only; consider MIPFocus=1 rerun
       INFEASIBLE                   -> computeIIS(); inspect conflicting rows/bounds
       INF_OR_UNBD                  -> set DualReductions=0, re-optimize to disambiguate
       UNBOUNDED                    -> missing constraint or sign error in the objective
7. validate the extracted solution with code that never touches the model

Construction idioms first — addVars index sets, tupledict.sum, generator-based addConstrs with named families, quicksum objectives — on a small transportation model:

"""Core gurobipy construction idioms on a tiny transportation model."""
import gurobipy as gp
from gurobipy import GRB

plants = ["p1", "p2"]
markets = ["m1", "m2", "m3"]
supply = {"p1": 30.0, "p2": 25.0}
demand = {"m1": 15.0, "m2": 20.0, "m3": 18.0}
cost = {
    ("p1", "m1"): 4.0, ("p1", "m2"): 6.0, ("p1", "m3"): 9.0,
    ("p2", "m1"): 5.0, ("p2", "m2"): 3.0, ("p2", "m3"): 7.0,
}

model = gp.Model("transport_idioms")
model.Params.OutputFlag = 0

# addVars over two index sets -> tupledict keyed by (plant, market) tuples.
ship = model.addVars(plants, markets, lb=0.0, name="ship")

# tupledict.sum with wildcards replaces hand-written quicksum for pure sums.
model.addConstrs((ship.sum(i, "*") <= supply[i] for i in plants), name="supply")
model.addConstrs((ship.sum("*", j) >= demand[j] for j in markets), name="demand")

model.setObjective(
    gp.quicksum(cost[i, j] * ship[i, j] for i in plants for j in markets),
    GRB.MINIMIZE,
)
model.optimize()

if model.Status == GRB.OPTIMAL:
    flows = {key: var.X for key, var in ship.items() if var.X > 1e-6}
    print(f"cost={model.ObjVal:.1f}")
    print(flows)
# Expected: cost=272.0 — p2 serves all of m2 (its comparative advantage), the
# remainder of p2 goes to m3, and p1 covers m1 plus the rest of m3.

The reusable harness wraps steps 5-7. It is solver-status-complete: it never reads .X without a usable incumbent, it disambiguates INF_OR_UNBD, and it returns a plain dataclass so downstream code (validators, result tables, plots) never holds gurobipy objects.

"""Reusable solve-and-extract harness with explicit status handling."""
from dataclasses import dataclass, field
from typing import Any

import gurobipy as gp
from gurobipy import GRB

STATUS_NAMES: dict[int, str] = {
    GRB.OPTIMAL: "OPTIMAL",
    GRB.INFEASIBLE: "INFEASIBLE",
    GRB.UNBOUNDED: "UNBOUNDED",
    GRB.INF_OR_UNBD: "INF_OR_UNBD",
    GRB.TIME_LIMIT: "TIME_LIMIT",
    GRB.INTERRUPTED: "INTERRUPTED",
}


@dataclass
class SolveResult:
    """Everything worth keeping from one solve, decoupled from gurobipy objects."""

    status_name: str
    objective: float | None
    best_bound: float | None
    mip_gap: float | None
    runtime_s: float
    solution: dict[str, dict[Any, float]] = field(default_factory=dict)

    @property
    def has_solution(self) -> bool:
        """True when objective and variable values are safe to use."""
        return self.objective is not None


def solve_and_extract(
    model: gp.Model,
    var_groups: dict[str, gp.tupledict],
    time_limit_s: float = 60.0,
    mip_gap: float = 1e-4,
    seed: int = 0,
) -> SolveResult:
    """Set parameters explicitly, optimize, extract values only if an incumbent exists."""
    model.Params.TimeLimit = time_limit_s
    model.Params.MIPGap = mip_gap
    model.Params.Seed = seed
    model.optimize()

    if model.Status == GRB.INF_OR_UNBD:
        # Presolve could not separate the two cases; rerun without dual reductions.
        model.Params.DualReductions = 0
        model.optimize()

    status_name = STATUS_NAMES.get(model.Status, f"STATUS_{model.Status}")
    usable = model.Status == GRB.OPTIMAL or (
        model.Status in (GRB.TIME_LIMIT, GRB.INTERRUPTED) and model.SolCount > 0
    )
    if not usable:
        return SolveResult(status_name, None, None, None, model.Runtime)

    is_mip = model.IsMIP == 1
    return SolveResult(
        status_name=status_name,
        objective=model.ObjVal,
        best_bound=model.ObjBound if is_mip else model.ObjVal,
        mip_gap=model.MIPGap if is_mip else 0.0,
        runtime_s=model.Runtime,
        solution={
            group: {key: var.X for key, var in tdict.items()}
            for group, tdict in var_groups.items()
        },
    )


# Tiny demo: a three-item 0-1 knapsack.
demo = gp.Model("harness_demo")
demo.Params.OutputFlag = 0
items = [0, 1, 2]
value = {0: 8.0, 1: 11.0, 2: 6.0}
weight = {0: 5.0, 1: 7.0, 2: 4.0}
pick = demo.addVars(items, vtype=GRB.BINARY, name="pick")
demo.addConstr(gp.quicksum(weight[i] * pick[i] for i in items) <= 10.0, name="budget")
demo.setObjective(gp.quicksum(value[i] * pick[i] for i in items), GRB.MAXIMIZE)

result = solve_and_extract(demo, {"pick": pick})
print(result.status_name, result.objective, result.solution["pick"])
# Expected: OPTIMAL 14.0 — items 0 and 2 chosen (weight 9 <= 10); item 1 alone gives only 11.

Worked Example 1: Production Mix — LP to MIP

The canonical walkthrough: start with the continuous production-mix LP, read its full dual story, then watch integrality and fixed charges change the answer. Products $j$ with unit profit $p_j$ compete for resources $r$ with capacity $b_r$, consuming $a_{rj}$ per unit:

$$ \max \sum_{j} p_j x_j \quad \text{s.t.} \quad \sum_{j} a_{rj}, x_j \le b_r ;; \forall r, \qquad x_j \ge 0. $$

The LP's dual values price the resources: $\pi_r$ is the marginal profit of one extra unit of resource $r$ (valid within the RHS sensitivity range), and a nonbasic product's reduced cost says how far its profit must rise before producing it becomes attractive.

"""Production-mix LP with full dual analysis: shadow prices, reduced costs, ranges."""
import gurobipy as gp
from gurobipy import GRB

PRODUCTS = ["chair", "table", "desk"]
RESOURCES = ["wood", "labor"]
PROFIT = {"chair": 45.0, "table": 80.0, "desk": 110.0}
CAPACITY = {"wood": 400.0, "labor": 450.0}
USAGE = {  # (resource, product) -> units consumed per unit produced
    ("wood", "chair"): 5.0, ("wood", "table"): 20.0, ("wood", "desk"): 25.0,
    ("labor", "chair"): 10.0, ("labor", "table"): 15.0, ("labor", "desk"): 20.0,
}


def build_production_lp() -> tuple[gp.Model, gp.tupledict]:
    """Continuous production mix: max p'x subject to Ax <= b, x >= 0."""
    model = gp.Model("production_mix_lp")
    model.Params.OutputFlag = 0
    x = model.addVars(PRODUCTS, lb=0.0, name="x")
    model.addConstrs(
        (
            gp.quicksum(USAGE[r, p] * x[p] for p in PRODUCTS) <= CAPACITY[r]
            for r in RESOURCES
        ),
        name="capacity",
    )
    model.setObjective(gp.quicksum(PROFIT[p] * x[p] for p in PRODUCTS), GRB.MAXIMIZE)
    return model, x


model, x = build_production_lp()
model.optimize()
assert model.Status == GRB.OPTIMAL

print(f"objective = {model.ObjVal:.2f}")
for p in PRODUCTS:
    print(f"  x[{p}] = {x[p].X:7.3f}   reduced cost = {x[p].RC:7.3f}")
for r in RESOURCES:
    row = model.getConstrByName(f"capacity[{r}]")
    print(
        f"  {r}: shadow price = {row.Pi:.4f}, slack = {row.Slack:.2f}, "
        f"RHS valid in [{row.SARHSLow:.1f}, {row.SARHSUp:.1f}]"
    )
# Expected: objective = 2258.33 with x = (chair 21.667, table 0, desk 11.667).
# Table's reduced cost is -4.167: its profit must rise above 84.17 to enter the mix.
# Shadow prices: wood 1.3333, labor 3.8333 — labor is the more valuable bottleneck.

Now the MIP step. Real production runs in integer lots, a production line incurs a fixed setup cost $f_j$ when used, and a used line must produce at least a minimum lot $\ell_j$. Introduce binaries $y_j \in {0,1}$ and link them to the quantities:

$$ \max \sum_j \left( p_j x_j - f_j y_j \right) \quad \text{s.t.} \quad \sum_j a_{rj} x_j \le b_r ;\forall r, \qquad \ell_j, y_j \le x_j \le u_j, y_j ;\forall j, \qquad x_j \in \mathbb{Z}_{\ge 0},; y_j \in {0,1}. $$

The linking coefficient $u_j$ is the big-M, and it must be tight: the best data-derived bound is $u_j = \min_r \lfloor b_r / a_{rj} \rfloor$, the most production any single resource allows. A lazy 1e6 here would weaken the relaxation and slow branching for no benefit.

"""Production mix as a MIP: integer lots, fixed setup charges, minimum lot sizes."""
import math

import gurobipy as gp
from gurobipy import GRB

PRODUCTS = ["chair", "table", "desk"]
RESOURCES = ["wood", "labor"]
PROFIT = {"chair": 45.0, "table": 80.0, "desk": 110.0}
CAPACITY = {"wood": 400.0, "labor": 450.0}
USAGE = {
    ("wood", "chair"): 5.0, ("wood", "table"): 20.0, ("wood", "desk"): 25.0,
    ("labor", "chair"): 10.0, ("labor", "table"): 15.0, ("labor", "desk"): 20.0,
}
SETUP_COST = {"chair": 200.0, "table": 150.0, "desk": 300.0}
MIN_LOT = {"chair": 5, "table": 5, "desk": 5}


def build_production_mip() -> tuple[gp.Model, gp.tupledict, gp.tupledict]:
    """Integer lots with fixed charges; big-M values derived from data, never guessed."""
    model = gp.Model("production_mix_mip")
    # Tightest valid upper bound on x[p]: the most any single resource allows.
    max_lot = {
        p: math.floor(min(CAPACITY[r] / USAGE[r, p] for r in RESOURCES))
        for p in PRODUCTS
    }
    x = model.addVars(
        PRODUCTS, vtype=GRB.INTEGER, lb=0,
        ub=[max_lot[p] for p in PRODUCTS], name="x",
    )
    y = model.addVars(PRODUCTS, vtype=GRB.BINARY, name="y")

    model.addConstrs(
        (
            gp.quicksum(USAGE[r, p] * x[p] for p in PRODUCTS) <= CAPACITY[r]
            for r in RESOURCES
        ),
        name="capacity",
    )
    model.addConstrs((x[p] <= max_lot[p] * y[p] for p in PRODUCTS), name="setup_link")
    model.addConstrs((x[p] >= MIN_LOT[p] * y[p] for p in PRODUCTS), name="min_lot")

    model.setObjective(
        gp.quicksum(PROFIT[p] * x[p] - SETUP_COST[p] * y[p] for p in PRODUCTS),
        GRB.MAXIMIZE,
    )
    return model, x, y


model, x, y = build_production_mip()
model.Params.OutputFlag = 0
model.Params.MIPGap = 1e-9  # prove optimality on this tiny instance
model.optimize()
assert model.Status == GRB.OPTIMAL

plan = {p: round(x[p].X) for p in PRODUCTS}
setups = {p: round(y[p].X) for p in PRODUCTS}
print(f"objective = {model.ObjVal:.1f}  plan = {plan}  setups = {setups}")
# Expected: objective = 1850.0 with plan chair=24, table=14, desk=0.
# The LP loved desks; the 300 setup charge flips the optimal mix to chairs + tables —
# integrality and fixed costs change structure, not just round the LP answer.

The qualitative lesson generalizes: never present a rounded LP solution as "near-optimal" when fixed charges or logical constraints are present. Here rounding the LP would keep desks open and lose 105 profit versus the true optimum, and on larger instances rounding is often infeasible outright.

Worked Example 2: Generalized Assignment Problem

GAP assigns each job $j$ to exactly one machine $i$, paying cost $c_{ij}$ and consuming capacity $a_{ij}$ of machine $i$'s budget $b_i$ (Ross & Soland (1975); Martello & Toth (1990), Knapsack Problems):

$$ \min \sum_{i}\sum_{j} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_{i} x_{ij} = 1 ;; \forall j, \qquad \sum_{j} a_{ij} x_{ij} \le b_i ;; \forall i, \qquad x_{ij} \in {0,1}. $$

The two constraint families have different characters — partitioning rows and knapsack rows — so they get one builder function each. This is exactly the structure that Lagrangian and column-generation methods later exploit; building it cleanly now pays off when the model becomes a subproblem.

"""Generalized assignment problem: min-cost job-to-machine assignment under capacities."""
from dataclasses import dataclass

import gurobipy as gp
from gurobipy import GRB


@dataclass(frozen=True)
class GapInstance:
    """GAP data: cost[i][j] and usage[i][j] for machine i and job j; capacity[i]."""

    cost: tuple[tuple[float, ...], ...]
    usage: tuple[tuple[float, ...], ...]
    capacity: tuple[float, ...]

    @property
    def machines(self) -> range:
        """Machine index set."""
        return range(len(self.capacity))

    @property
    def jobs(self) -> range:
        """Job index set."""
        return range(len(self.cost[0]))


def add_assignment_constraints(model: gp.Model, x: gp.tupledict, inst: GapInstance) -> None:
    """Each job goes to exactly one machine (set-partitioning rows)."""
    model.addConstrs(
        (gp.quicksum(x[i, j] for i in inst.machines) == 1 for j in inst.jobs),
        name="assign",
    )


def add_capacity_constraints(model: gp.Model, x: gp.tupledict, inst: GapInstance) -> None:
    """Per-machine resource budget (knapsack rows)."""
    model.addConstrs(
        (
            gp.quicksum(inst.usage[i][j] * x[i, j] for j in inst.jobs)
            <= inst.capacity[i]
            for i in inst.machines
        ),
        name="capacity",
    )


def build_gap_model(inst: GapInstance) -> tuple[gp.Model, gp.tupledict]:
    """Assemble the GAP MIP: one builder function per constraint family."""
    model = gp.Model("gap")
    x = model.addVars(inst.machines, inst.jobs, vtype=GRB.BINARY, name="x")
    add_assignment_constraints(model, x, inst)
    add_capacity_constraints(model, x, inst)
    model.setObjective(
        gp.quicksum(
            inst.cost[i][j] * x[i, j] for i in inst.machines for j in inst.jobs
        ),
        GRB.MINIMIZE,
    )
    return model, x


def solve_gap(inst: GapInstance, time_limit_s: float = 30.0) -> dict[int, int] | None:
    """Solve and return job -> machine, or None when no incumbent exists."""
    model, x = build_gap_model(inst)
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit_s
    model.Params.MIPGap = 1e-9
    model.optimize()
    usable = model.Status == GRB.OPTIMAL or (
        model.Status == GRB.TIME_LIMIT and model.SolCount > 0
    )
    if not usable:
        return None
    print(f"objective = {model.ObjVal:.1f}  gap = {model.MIPGap:.2%}")
    return {j: i for i in inst.machines for j in inst.jobs if x[i, j].X > 0.5}


INST = GapInstance(
    cost=((6.0, 9.0, 4.0, 8.0), (5.0, 3.0, 8.0, 6.0)),
    usage=((3.0, 4.0, 2.0, 5.0), (4.0, 2.0, 5.0, 3.0)),
    capacity=(8.0, 8.0),
)
print(solve_gap(INST))
# Expected: objective = 19.0, assignment {0: 0, 1: 1, 2: 0, 3: 1}.
# Cheapest row-wise assignment costs 18 but overloads machine 1 (load 9 > 8);
# moving job 0 to machine 0 is the cheapest repair (+1).

Every model deserves an independent validator: a function that recomputes feasibility and the objective from raw data and a plain solution, importing nothing from the solver. It catches indexing bugs, wrong-direction inequalities, and silent unit errors that the solver cannot see (it faithfully optimizes whatever wrong model you gave it).

"""Independent feasibility and objective check for a GAP solution (no solver involved)."""


def validate_gap_solution(
    assignment: dict[int, int],
    cost: list[list[float]],
    usage: list[list[float]],
    capacity: list[float],
    tol: float = 1e-6,
) -> tuple[float, list[str]]:
    """Recompute objective and collect violations; an empty list means feasible."""
    n_machines, n_jobs = len(capacity), len(cost[0])
    violations: list[str] = []
    if sorted(assignment) != list(range(n_jobs)):
        violations.append("not every job is assigned exactly once")
    total = 0.0
    load = [0.0] * n_machines
    for job, machine in assignment.items():
        total += cost[machine][job]
        load[machine] += usage[machine][job]
    for i in range(n_machines):
        if load[i] > capacity[i] + tol:
            violations.append(f"machine {i} overloaded: {load[i]} > {capacity[i]}")
    return total, violations


obj, errors = validate_gap_solution(
    {0: 0, 1: 1, 2: 0, 3: 1},
    cost=[[6.0, 9.0, 4.0, 8.0], [5.0, 3.0, 8.0, 6.0]],
    usage=[[3.0, 4.0, 2.0, 5.0], [4.0, 2.0, 5.0, 3.0]],
    capacity=[8.0, 8.0],
)
print(obj, errors)
# Expected: 19.0 [] — matches the solver's objective with zero violations.

Cross-check the validator's objective against model.ObjVal on every run; a mismatch larger than 1e-6 * |ObjVal| means the model and your intent have diverged.

Advanced Techniques

Matrix API for large, regular structure

When the model is naturally expressed with arrays — knapsacks, assignment blocks, flow balances over dense index grids — addMVar plus @ builds it one or two orders of magnitude faster than per-variable loops, because construction happens in C on whole arrays.

"""Matrix API: build a model from numpy arrays with addMVar."""
import numpy as np

import gurobipy as gp
from gurobipy import GRB

value = np.array([8.0, 11.0, 6.0, 4.0])
weight = np.array([5.0, 7.0, 4.0, 3.0])
budget = 10.0

model = gp.Model("knapsack_matrix")
model.Params.OutputFlag = 0
x = model.addMVar(value.size, vtype=GRB.BINARY, name="x")
model.addConstr(weight @ x <= budget, name="budget")
model.setObjective(value @ x, GRB.MAXIMIZE)
model.optimize()

assert model.Status == GRB.OPTIMAL
chosen = np.flatnonzero(x.X > 0.5)
print(f"objective = {model.ObjVal:.1f}  chosen = {chosen.tolist()}")
# Expected: objective = 15.0, chosen = [1, 3] (weights 7 + 3 = 10, exactly the budget).

For block-structured constraints pass a scipy.sparse matrix to model.addMConstr or use A @ x <= b directly with sparse A. Mixing styles is fine: use addMVar for the bulk arrays and addConstrs generators for the irregular side constraints.

Parameter strategy by situation

Set parameters explicitly and record them; defaults change between Gurobi versions, and an unrecorded parameter is an unreproducible experiment.

"""Named parameter profiles: explicit, recorded, reusable."""
import gurobipy as gp

PROFILES: dict[str, dict[str, float | int]] = {
    "prove_optimality": {"MIPGap": 1e-9, "MIPFocus": 2},
    "good_solution_fast": {"MIPGap": 0.01, "MIPFocus": 1, "Heuristics": 0.2, "TimeLimit": 60},
    "hard_root_lp": {"Method": 2, "Crossover": 0},  # barrier, skip crossover
    "reproducible_run": {"Seed": 42, "Threads": 1},
    "numerically_fragile": {"NumericFocus": 2, "IntFeasTol": 1e-7},
}


def apply_profile(model: gp.Model, profile: str) -> None:
    """Apply a named parameter bundle to a model."""
    for param, value in PROFILES[profile].items():
        model.setParam(param, value)


demo = gp.Model("profile_demo")
demo.Params.OutputFlag = 0
apply_profile(demo, "good_solution_fast")
print(demo.Params.MIPGap, demo.Params.TimeLimit)
# Expected: 0.01 60.0

Situation	Parameters that matter
Need any good solution quickly	MIPFocus=1, Heuristics=0.2, modest MIPGap (1-5%)
Bound moves too slowly	MIPFocus=3, Cuts=2
Proving optimality at the end	MIPFocus=2, tight MIPGap
Root LP dominates runtime	Method=2 (barrier), Crossover=0 if a basis is not needed
Benchmarking / experiments	Seed, Threads=1, fixed TimeLimit — see the run protocol in your experiment harness
Memory pressure on big trees	NodefileStart=0.5 to spill nodes to disk

Diagnosing infeasibility with an IIS

An Irreducible Inconsistent Subsystem is a minimal set of constraints and bounds that is infeasible together — removing any one member makes it feasible. It turns "model infeasible" into a short, readable list of suspects. Named constraints make the output meaningful.

"""Infeasibility diagnosis: compute and print an IIS."""
import gurobipy as gp
from gurobipy import GRB

model = gp.Model("infeasible_demo")
model.Params.OutputFlag = 0
x = model.addVar(lb=0.0, name="x")
y = model.addVar(lb=0.0, name="y")
model.addConstr(x + y >= 10.0, name="demand")
model.addConstr(x <= 3.0, name="cap_x")
model.addConstr(y <= 4.0, name="cap_y")
model.optimize()

if model.Status == GRB.INFEASIBLE:
    model.computeIIS()
    members = [c.ConstrName for c in model.getConstrs() if c.IISConstr]
    print("IIS:", members)
# Expected: IIS: ['demand', 'cap_x', 'cap_y'] — max supply 3 + 4 = 7 < 10 demanded;
# dropping any one of the three rows restores feasibility.

For "which constraints should I relax and by how much," model.feasRelax finds a minimal weighted relaxation. Both tools, plus callbacks and the solution pool, are covered in depth in gurobi-advanced-features.

Numerical hygiene

Keep all matrix coefficients within roughly six orders of magnitude (ideally 1e-3 to 1e6), rescale units rather than fight tolerances, and derive every big-M from data as in the production MIP. Symptoms of trouble: warnings about large coefficient ranges, "solutions" that your independent validator rejects, and binaries at 0.5 after "optimal" termination. Klotz & Newman (2013) is the standard practical reference. IntFeasTol (default 1e-5) means a binary can legitimately return 0.99999 — always extract with > 0.5, and tighten the tolerance rather than post-round when downstream logic is sensitive.

Warm starts and re-solve loops

Set var.Start = value before optimize() to inject a known feasible solution; Gurobi repairs partial starts. In rolling-horizon or decomposition loops, modify RHS/Obj/bounds in place instead of rebuilding — the retained basis or incumbent typically cuts re-solve time substantially. Feed heuristic solutions in as MIP starts whenever you have them; a good incumbent at the root prunes the tree from the first node.

Practical Challenges

Model construction takes longer than the solve. Python-side expression building is the usual culprit: sum(...) over LinExpr objects is quadratic, one addConstr call per row inside nested loops pays call overhead millions of times. Use gp.quicksum, generator-based addConstrs, and the matrix API for the dense blocks. Profile build time separately from model.Runtime — they are different problems with different fixes.

AttributeError (unable to retrieve attribute X). You read a solution that does not exist: before optimize(), after INFEASIBLE, or after TIME_LIMIT with SolCount == 0. Gate every extraction on status as in the harness; never wrap the read in a bare try/except that hides which case occurred.

Binaries come back as 0.9999999 or 1e-9. Integrality is enforced only to IntFeasTol. Test binaries with > 0.5, recover integer counts with round(var.X), and never use equality comparisons on solution values. If downstream feasibility is sensitive to the residue, tighten IntFeasTol to 1e-7 and re-solve.

The model is infeasible but the data "looks fine." Compute an IIS and read the named members. The most common root causes, in order: an inequality written in the wrong direction, a variable's default lb=0 clashing with a constraint that needs it negative, unit mismatches between two data sources, and over-tight bounds copied from another instance. model.write("model.lp") and reading the rows around the IIS members settles it quickly.

Big-M chosen as 1e6 "to be safe." Oversized big-M weakens the LP relaxation (the solver branches far more) and creates numerical traps: with FeasibilityTol = 1e-6, a big-M of 1e6 lets y = 1e-12 "switch off" a constraint while the incumbent quietly violates your intent. Derive each big-M from capacities and bounds, per entity, as in the production MIP; if no finite bound exists, reconsider the formulation (indicator constraints are the fallback — see linearization-techniques).

The solver says optimal but the answer is wrong in the real world. The solver proved optimality of the model you wrote, not the problem you meant. Run the independent validator, recompute the objective from raw data, and compare term by term. Divergence localizes the bug to a constraint family or a cost term; this is why builders are separate functions you can test one at a time.

The MIP gap stalls for hours. This is usually formulation weakness or symmetry, not a parameter issue: identical machines, interchangeable bins, or loose linking constraints leave the solver enumerating permutations of the same solution. Tighten the relaxation, add symmetry-breaking, or reformulate — the diagnosis-and-repair playbook is integer-programming-techniques. Parameter nudges (MIPFocus, Cuts) buy at most a constant factor.

Every data change triggers a full model rebuild. For re-solve loops, mutate the model: constr.RHS = new_value, var.Obj = new_cost, var.UB = 0.0 to forbid, model.chgCoeff for matrix entries. Rebuild only when the sparsity pattern itself changes. This preserves warm-start information and often turns minutes per iteration into seconds.

Tools & Libraries

Library	When to use	Note
gurobipy	The subject of this skill	Commercial license required; free academic licenses; pip install includes a size-limited trial
numpy	Data prep and the matrix API	addMVar + @ is the fast path for dense regular blocks
scipy.sparse	Constraint matrices for addMConstr	Build COO/CSR once; avoid dense matrices above ~1e4 columns
pandas	Solution and experiment tables	Convert SolveResult objects to tidy one-row-per-run tables
PuLP / Pyomo / python-mip	Solver-agnostic modeling layers	Trade API expressiveness for portability — see open-source-solvers
HiGHS / SCIP / CBC	No Gurobi license available	Selection and migration guidance in open-source-solvers
OR-Tools CP-SAT	Disjunctive scheduling, heavy logical structure	Different modeling paradigm; triage in problem-formulation

Output Format

A complete MILP modeling deliverable contains, in order:

1. Formulation. Sets, parameters, decision variables (with types and bounds), constraints, and objective in mathematical notation — before any code. Every constraint family in the code must point back to a numbered line here.

2. Model summary table. What was actually built:

   Component	Count	Detail
   Binary variables	m·n	x[i,j] assignment
   Integer variables	n	x[p], ub from capacity
   Constraint family assign	n	partitioning rows
   Constraint family capacity	m	knapsack rows

3. Solver configuration. Gurobi version, every parameter changed from default, seed, thread count, hardware if runtimes are reported.

4. Solution-quality report. Status, incumbent objective, best bound, relative gap, runtime, node count. For TIME_LIMIT results, state the gap explicitly — "best found 1850 (gap 2.3% after 600 s)" — never present an incumbent as optimal.

5. Decision tables. Nonzero decisions only (var.X > 0.5 for binaries, > 1e-6 for continuous), in domain terms ("job 3 on machine 1"), not raw variable dumps.

6. Validation statement. Output of the independent checker: recomputed objective, maximum constraint violation, and confirmation that it matches ObjVal within tolerance.

7. Artifacts when useful. model.write("model.lp") for review on small instances, model.write("model.mps") plus the data seed for reproducible bug reports, solution JSON/CSV for downstream consumers.

Questions to Ask

• Is the mathematical formulation written down, or are we starting from a word problem?
• How large is a realistic instance — how many variables and constraints as a function of the data?
• Which constraints are hard, and which are preferences that can be violated at a price?
• Is a Gurobi license available (academic or commercial), or should this target an open-source solver?
• What is the time budget per solve, and is proven optimality required or is a small gap acceptable?
• Will the model be solved once, or repeatedly with changing data (rolling horizon, decomposition loop)?
• Where do the coefficients come from, and are units consistent across data sources?
• Are there quantities currently modeled as integer that could safely be continuous?
• What should the deliverable look like — a solution table, a reusable module, or an experiment harness?

Related Skills

• linearization-techniques — when the formulation contains products of variables, absolute values, min/max terms, piecewise-linear costs, or logical implications that must become linear before this skill's machinery applies.
• integer-programming-techniques — when the model is correct but slow: relaxation strength, symmetry breaking, big-M vs indicator trade-offs, and MIP gap interpretation.
• gurobi-advanced-features — when you need callbacks (lazy constraints, user cuts, heuristic injection), IIS workflows, the solution pool, multi-objective API, or the parameter tuning tool.
• problem-formulation — when decisions, objective, and constraints are not yet pinned down, or when it is unclear whether MILP is the right model class at all.
• open-source-solvers — when no Gurobi license is available and the model should be built for or migrated to HiGHS, SCIP, CBC, or OR-Tools.