Mixed-Integer Linear Programming (MILP)#
A tutorial for first-year PhD students on formulating and solving mixed-integer linear programs.
Learning Objectives#
By the end of this notebook you will be able to:
Explain why adding integrality constraints to a linear program makes it NP-hard
Describe the Branch-and-Bound algorithm and its key concepts (incumbent, bound, gap, pruning)
Formulate three classic MILPs: 0-1 knapsack, facility location, and job-shop scheduling
Compare big-M constraints with indicator constraints and understand their impact on solver performance
Apply symmetry-breaking constraints to reduce the search space
Use SOS1 constraints for piecewise-linear modeling
Prerequisites: Basic linear programming, familiarity with Python and NumPy.
References: This tutorial draws on the textbooks by [Wolsey, 1998] and [Nemhauser and Wolsey, 1988].
import discopt.modeling as dm
import numpy as np
1. From LP to MILP#
A linear program (LP) optimizes a linear objective subject to linear constraints over continuous variables:
LPs are solvable in polynomial time (via the simplex or interior-point methods). But many real-world decisions are discrete: build a factory or not, assign a job to a machine, select a subset of items.
A mixed-integer linear program (MILP) adds integrality constraints on some or all variables:
This seemingly small change has profound consequences:
The feasible region is no longer convex (it is a finite set of points)
The problem becomes NP-hard in general
Solving requires combinatorial search, not just linear algebra
The dominant algorithm for solving MILPs is Branch and Bound [Land and Doig, 1960], first proposed by Land and Doig in 1960 and refined by Dakin [Dakin, 1965].
2. Branch-and-Bound Overview#
Branch and Bound (B&B) systematically explores the space of integer variable assignments by solving a tree of LP relaxations.
Key concepts:
Term |
Definition |
|---|---|
LP relaxation |
Drop integrality requirements; solve the continuous LP |
Lower bound |
The LP relaxation objective (for minimization) |
Upper bound (incumbent) |
Best feasible integer solution found so far |
Gap |
$(\text{UB} - \text{LB}) / |
Branching |
Pick a fractional integer variable \(y_i = 0.7\), create two children: \(y_i \leq 0\) and \(y_i \geq 1\) |
Pruning |
Discard a node if its LB \(\geq\) UB (cannot improve on incumbent) |
Fathoming |
A node is fathomed if it is pruned, infeasible, or integer-feasible |
B&B Tree (ASCII illustration)#
[Root: LP relaxation]
obj = 3.7, y1 = 0.6
/ \
y1 <= 0 y1 >= 1
/ \
[Node 1] [Node 2]
obj = 4.2 obj = 3.9, y2 = 0.4
integer feasible! / \
UB = 4.2 y2 <= 0 y2 >= 1
/ \
[Node 3] [Node 4]
obj = 4.5 obj = 4.0
PRUNED (4.5 > 4.2) integer feasible!
UB = 4.0 (new best!)
The algorithm terminates when all nodes are fathomed. At that point, the gap is zero and the incumbent is provably optimal.
In discopt, the B&B tree is managed by a Rust backend for efficiency, and the per-node LP relaxations are solved by the in-house Rust simplex. A model that classifies as a pure MILP skips the spatial machinery entirely and is routed to HiGHS, with discopt re-verifying the certificate.
3. Example 1: 0-1 Knapsack Problem#
The classic combinatorial optimization problem: given \(n\) items, each with a weight \(w_i\) and a value \(v_i\), select a subset that maximizes total value without exceeding a weight capacity \(W\).
This is one of Karp’s 21 NP-complete problems.
# Problem data: 10 items
weights = np.array([10, 20, 30, 40, 50, 15, 25, 35, 45, 5])
values = np.array([60, 100, 120, 140, 150, 50, 80, 110, 130, 20])
capacity = 50
n_items = len(weights)
print(f"Items: {n_items}")
print(f"Capacity: {capacity}")
print(f"Value-to-weight ratios: {values / weights}")
Items: 10
Capacity: 50
Value-to-weight ratios: [6. 5. 4. 3.5 3. 3.33333333
3.2 3.14285714 2.88888889 4. ]
m = dm.Model("knapsack")
# Binary decision: select item i or not
x = m.binary("x", shape=(n_items,))
# Maximize total value
m.maximize(dm.sum(lambda i: values[i] * x[i], over=range(n_items)))
# Weight capacity constraint
m.subject_to(dm.sum(lambda i: weights[i] * x[i], over=range(n_items)) <= capacity, name="capacity")
print(m.summary())
Model: knapsack
Variables: 10 (0 continuous, 10 integer/binary)
Constraints: 1
Objective: maximize Σ[10 terms]
Parameters: 0
result = m.solve()
print(f"Status: {result.status}")
print(f"Objective: ${result.objective:.2f}")
print(f"Gap: {result.gap}")
print(f"Nodes: {result.node_count}")
print(f"Time: {result.wall_time:.3f}s")
print()
# Display selected items
selected = result.value(x)
print("Selected items:")
total_weight = 0
total_value = 0
for i in range(n_items):
if selected[i] > 0.5: # binary threshold
print(f" Item {i}: weight={weights[i]}, value={values[i]}")
total_weight += weights[i]
total_value += values[i]
print(f"\nTotal weight: {total_weight} / {capacity}")
print(f"Total value: {total_value}")
Status: optimal
Objective: $230.00
Gap: 0.0
Nodes: 1
Time: 0.040s
Selected items:
Item 0: weight=10, value=60
Item 1: weight=20, value=100
Item 5: weight=15, value=50
Item 9: weight=5, value=20
Total weight: 50 / 50
Total value: 230
Interpreting the results#
Gap = 0 means the solution is provably optimal. The B&B algorithm has explored enough nodes to certify that no better integer solution exists.
Node count tells you how hard the problem was. A small node count means the LP relaxation was tight (close to the integer optimum) and many branches were pruned early.
The knapsack LP relaxation is famously tight: it typically selects all items by decreasing value-to-weight ratio and splits only the last item fractionally.
4. Example 2: Uncapacitated Facility Location#
A company must decide which of \(m\) candidate facilities to open in order to serve \(n\) customers at minimum total cost (fixed opening costs plus shipping costs).
subject to: $\(\sum_{i=1}^{m} x_{ij} \geq 1 \quad \forall j \quad \text{(each customer served)}\)\( \)\(x_{ij} \leq y_i \quad \forall i, j \quad \text{(only open facilities can serve)}\)\( \)\(y_i \in \{0, 1\}, \; x_{ij} \in [0, 1]\)$
Here \(y_i = 1\) if facility \(i\) is opened, and \(x_{ij}\) is the fraction of customer \(j\)’s demand served by facility \(i\).
# Problem data
n_facilities = 6
n_customers = 15
# Fixed costs to open each facility
fixed_costs = np.array([100, 120, 90, 110, 130, 95])
# Customer demands
demands = np.array([10, 15, 12, 8, 20, 18, 14, 9, 11, 16, 13, 7, 19, 17, 6])
# Shipping cost per unit from facility i to customer j
np.random.seed(42)
shipping_costs = np.random.randint(1, 20, size=(n_facilities, n_customers)).astype(float)
print(f"Facilities: {n_facilities}")
print(f"Customers: {n_customers}")
print(f"Fixed costs: {fixed_costs}")
print(f"Total demand: {demands.sum()}")
Facilities: 6
Customers: 15
Fixed costs: [100 120 90 110 130 95]
Total demand: 195
# Formulation with big-M linking constraints: x[i,j] <= y[i]
m = dm.Model("facility_location_bigM")
# Binary: open facility i?
y = m.binary("open", shape=(n_facilities,))
# Continuous: fraction of customer j's demand served by facility i
x = m.continuous("alloc", shape=(n_facilities, n_customers), lb=0, ub=1)
# Objective: fixed costs + shipping costs
m.minimize(
dm.sum(lambda i: fixed_costs[i] * y[i], over=range(n_facilities))
+ dm.sum(
lambda i: dm.sum(
lambda j: shipping_costs[i, j] * demands[j] * x[i, j], over=range(n_customers)
),
over=range(n_facilities),
)
)
# Each customer must be fully served
for j in range(n_customers):
m.subject_to(dm.sum(lambda i: x[i, j], over=range(n_facilities)) >= 1, name=f"demand_{j}")
# Big-M linking: can only allocate from open facilities
for i in range(n_facilities):
for j in range(n_customers):
m.subject_to(x[i, j] <= y[i], name=f"link_{i}_{j}")
result_bigm = m.solve()
print(f"Status: {result_bigm.status}")
print(f"Objective: {result_bigm.objective:.2f}")
print(f"Gap: {result_bigm.gap}")
print(f"Nodes: {result_bigm.node_count}")
print(f"Time: {result_bigm.wall_time:.3f}s")
print()
# Which facilities are open?
open_facilities = result_bigm.value(y)
for i in range(n_facilities):
status = "OPEN" if open_facilities[i] > 0.5 else "closed"
print(f" Facility {i} (cost={fixed_costs[i]}): {status}")
Status: optimal
Objective: 1128.00
Gap: 0.0
Nodes: 1
Time: 0.017s
Facility 0 (cost=100): OPEN
Facility 1 (cost=120): closed
Facility 2 (cost=90): OPEN
Facility 3 (cost=110): OPEN
Facility 4 (cost=130): closed
Facility 5 (cost=95): OPEN
Big-M vs Indicator Constraints#
The constraint \(x_{ij} \leq y_i\) is already the tightest possible big-M linking for this formulation (since \(x_{ij} \leq 1\) and \(M = 1\)). In other problems, choosing a good value of \(M\) is critical:
\(M\) too large: The LP relaxation becomes very weak, leading to a large gap and many B&B nodes.
\(M\) too small: The formulation may cut off feasible solutions.
Indicator constraints avoid the big-M entirely. Instead of \(x_{ij} \leq M \cdot y_i\), we write:
“If \(y_i = 1\), then the allocation constraints for facility \(i\) must hold.”
discopt supports indicator constraints via m.if_then(indicator, then_constraints),
which enforces then_constraints only when the binary indicator is 1. Note the
direction: here we need “if facility \(i\) is closed, its allocations are zero”,
so the indicator is the complement \(1 - y_i\), not \(y_i\) itself.
# Reformulation with genuine indicator constraints (m.if_then)
m2 = dm.Model("facility_location_indicator")
y2 = m2.binary("open", shape=(n_facilities,))
x2 = m2.continuous("alloc", shape=(n_facilities, n_customers), lb=0, ub=1)
# Same objective
m2.minimize(
dm.sum(lambda i: fixed_costs[i] * y2[i], over=range(n_facilities))
+ dm.sum(
lambda i: dm.sum(
lambda j: shipping_costs[i, j] * demands[j] * x2[i, j], over=range(n_customers)
),
over=range(n_facilities),
)
)
# Demand constraints (same as before)
for j in range(n_customers):
m2.subject_to(dm.sum(lambda i: x2[i, j], over=range(n_facilities)) >= 1, name=f"demand_{j}")
# Indicator constraints. m.if_then fires when its indicator is 1, and the logic we
# need is "if facility i is CLOSED, all of its allocations are zero" -- so the
# indicator is the complement of y2[i], introduced as its own binary.
closed = m2.binary("closed", shape=(n_facilities,))
for i in range(n_facilities):
m2.subject_to(closed[i] == 1 - y2[i], name=f"complement_{i}")
m2.if_then(
closed[i],
[x2[i, j] <= 0 for j in range(n_customers)],
name=f"shut_{i}",
)
result_indicator = m2.solve()
print(f"Status: {result_indicator.status}")
print(f"Objective: {result_indicator.objective:.2f}")
print(f"Gap: {result_indicator.gap}")
print(f"Nodes: {result_indicator.node_count}")
print(f"Time: {result_indicator.wall_time:.3f}s")
print()
print(f"Big-M objective: {result_bigm.objective:.2f}")
print(f"Indicator objective: {result_indicator.objective:.2f}")
print(f"Big-M nodes: {result_bigm.node_count}")
print(f"Indicator nodes: {result_indicator.node_count}")
Status: optimal
Objective: 1128.00
Gap: 0.0
Nodes: 1
Time: 0.024s
Big-M objective: 1128.00
Indicator objective: 1128.00
Big-M nodes: 1
Indicator nodes: 1
Discussion: Formulation Strength#
Both models above return the same optimum and solve in the same number of nodes. That is the expected outcome here, and it is worth stating plainly rather than glossing over: as noted above, \(x_{ij} \le y_i\) is already the tightest big-M for this problem (\(M = 1\)), so there is no looseness left for an indicator constraint to recover. Indicator constraints earn their keep when the natural big-M is much larger than the variable it multiplies — a start time bounded by the sum of all processing times, say, as in the job-shop model below.
A stronger formulation is one whose LP relaxation is closer to the convex hull of integer-feasible points. Key insights:
Tighter LP relaxation \(\Rightarrow\) better lower bounds \(\Rightarrow\) more pruning \(\Rightarrow\) fewer nodes
The “ideal” formulation would have its LP relaxation equal to the convex hull, but computing the convex hull is generally as hard as solving the MILP itself
In practice, good formulations strike a balance between tightness and the number of constraints/variables
This is a central theme in integer programming: the formulation matters as much as the algorithm [Wolsey, 1998].
5. Example 3: Job-Shop Scheduling#
Schedule \(n\) jobs on \(K\) machines to minimize the makespan (total completion time). Each job \(i\) has a processing time \(p_{ik}\) on machine \(k\). On each machine, jobs cannot overlap.
Variables:
\(s_{ik} \geq 0\): start time of job \(i\) on machine \(k\)
\(z_{ijk} \in \{0, 1\}\): 1 if job \(i\) is processed before job \(j\) on machine \(k\)
\(C_{\max}\): makespan
Disjunctive constraints (big-M): $\(s_{ik} + p_{ik} \leq s_{jk} + M(1 - z_{ijk}) \quad \forall i < j, \; \forall k\)\( \)\(s_{jk} + p_{jk} \leq s_{ik} + M \cdot z_{ijk} \quad \forall i < j, \; \forall k\)$
These ensure that for each pair of jobs on the same machine, one must finish before the other starts.
# Problem data: 5 jobs, 3 machines
n_jobs = 5
n_machines = 3
# Processing times: proc[i, k] = time for job i on machine k
proc = np.array(
[
[3, 2, 4], # Job 0
[2, 3, 1], # Job 1
[4, 1, 3], # Job 2
[1, 4, 2], # Job 3
[3, 2, 3], # Job 4
]
)
# Big-M: upper bound on any start time
M = float(proc.sum()) # sum of all processing times is a valid upper bound
print(f"Jobs: {n_jobs}, Machines: {n_machines}")
print(f"Processing times:\n{proc}")
print(f"Big-M value: {M}")
Jobs: 5, Machines: 3
Processing times:
[[3 2 4]
[2 3 1]
[4 1 3]
[1 4 2]
[3 2 3]]
Big-M value: 38.0
m = dm.Model("jobshop")
# Start times for each job on each machine
s = m.continuous("start", shape=(n_jobs, n_machines), lb=0, ub=M)
# Makespan variable
C = m.continuous("makespan", lb=0, ub=M)
# Ordering binaries: z[i,j,k] = 1 if job i before job j on machine k (i < j)
# We flatten to a 1D array for simplicity
n_pairs = n_jobs * (n_jobs - 1) // 2
z = m.binary("order", shape=(n_pairs * n_machines,))
# Minimize makespan
m.minimize(C)
# Makespan >= completion time of every job on every machine
for i in range(n_jobs):
for k in range(n_machines):
m.subject_to(C >= s[i, k] + proc[i, k], name=f"makespan_{i}_{k}")
# Disjunctive constraints: for each pair (i, j) with i < j on each machine k
pair_idx = 0
for i in range(n_jobs):
for j in range(i + 1, n_jobs):
for k in range(n_machines):
idx = pair_idx * n_machines + k
# If z = 1: job i before job j
m.subject_to(
s[i, k] + proc[i, k] <= s[j, k] + M * (1 - z[idx]), name=f"disj_fwd_{i}_{j}_{k}"
)
# If z = 0: job j before job i
m.subject_to(s[j, k] + proc[j, k] <= s[i, k] + M * z[idx], name=f"disj_bwd_{i}_{j}_{k}")
pair_idx += 1
print(
f"Variables: {n_jobs * n_machines} start times + 1 makespan + {n_pairs * n_machines} binaries"
)
print(f"Constraints: {n_jobs * n_machines} makespan + {2 * n_pairs * n_machines} disjunctive")
Variables: 15 start times + 1 makespan + 30 binaries
Constraints: 15 makespan + 60 disjunctive
result = m.solve(time_limit=60)
print(f"Status: {result.status}")
obj_str = f"{result.objective:.2f}" if result.objective is not None else "(no incumbent)"
print(f"Makespan: {obj_str}")
print(f"Gap: {result.gap}")
print(f"Nodes: {result.node_count}")
print(f"Time: {result.wall_time:.3f}s")
Status: optimal
Makespan: 13.00
Gap: 0.0
Nodes: 984
Time: 0.437s
# Text-based Gantt chart
starts = result.value(s)
makespan_val = result.objective
print("\nGantt Chart (each '=' represents 1 time unit):")
print("=" * 60)
for k in range(n_machines):
# Build a timeline for this machine
timeline_len = int(np.ceil(makespan_val)) + 1
timeline = ["."] * timeline_len
for i in range(n_jobs):
start_t = int(round(starts[i, k]))
for t in range(start_t, start_t + proc[i, k]):
if t < timeline_len:
timeline[t] = str(i)
bar = "".join(timeline)
print(f"Machine {k}: |{bar}|")
print("=" * 60)
print("Legend: digits = job index, '.' = idle")
print(f"Makespan = {makespan_val:.1f}")
Gantt Chart (each '=' represents 1 time unit):
============================================================
Machine 0: |4443222211000.|
Machine 1: |4433332.11100.|
Machine 2: |4443322210000.|
============================================================
Legend: digits = job index, '.' = idle
Makespan = 13.0
6. Symmetry Breaking (and how to get it wrong)#
Many MILPs have symmetric solutions: permuting interchangeable objects yields an equally good solution. B&B wastes effort exploring these equivalent branches, so we add symmetry-breaking constraints that impose an arbitrary ordering — for example \(y_0 \ge y_1 \ge \cdots\) (“if we open one, prefer the lower index”).
The word doing all the work in that sentence is interchangeable. A
symmetry-breaking constraint is only valid when permuting the objects genuinely
maps feasible solutions to equally-costly feasible solutions. If the objects differ
in any coefficient, the ordering constraint is not a symmetry break at all — it
is an invalid cut that can remove the true optimum, and the solver will happily
certify the wrong answer as optimal, because you told it that constraint was
part of the problem.
The next cell demonstrates both cases side by side.
n_fac = 4
n_cust = 10
identical_cost = 100.0
rng = np.random.RandomState(123)
# Two cost structures:
# - "distinct": facilities share a fixed cost but have DIFFERENT shipping costs,
# so they are NOT interchangeable.
# - "identical": every facility has the SAME shipping-cost row, so any permutation
# of the facilities is a genuine symmetry.
ship_distinct = rng.randint(2, 15, size=(n_fac, n_cust)).astype(float)
dem = rng.randint(5, 20, size=(n_cust,)).astype(float)
ship_identical = np.tile(ship_distinct[0], (n_fac, 1))
def build_facility_model(name, ship, add_symmetry_breaking=False):
"""Build a facility location model, optionally with an ordering constraint."""
mdl = dm.Model(name)
y = mdl.binary("open", shape=(n_fac,))
x = mdl.continuous("alloc", shape=(n_fac, n_cust), lb=0, ub=1)
mdl.minimize(
dm.sum(lambda i: identical_cost * y[i], over=range(n_fac))
+ dm.sum(
lambda i: dm.sum(lambda j: ship[i, j] * dem[j] * x[i, j], over=range(n_cust)),
over=range(n_fac),
)
)
for j in range(n_cust):
mdl.subject_to(dm.sum(lambda i: x[i, j], over=range(n_fac)) >= 1, name=f"demand_{j}")
for i in range(n_fac):
for j in range(n_cust):
mdl.subject_to(x[i, j] <= y[i], name=f"link_{i}_{j}")
if add_symmetry_breaking:
# Ordering constraints: y[0] >= y[1] >= y[2] >= y[3]
for i in range(n_fac - 1):
mdl.subject_to(y[i] >= y[i + 1], name=f"sym_break_{i}")
return mdl
print("A) Facilities that are NOT interchangeable (shipping costs differ)")
a1 = build_facility_model("distinct_no_sym", ship_distinct, False).solve()
a2 = build_facility_model("distinct_sym", ship_distinct, True).solve()
print(f" without ordering constraint: obj={a1.objective:.2f} nodes={a1.node_count}")
print(f" with ordering constraint: obj={a2.objective:.2f} nodes={a2.node_count}")
print(f" --> the ordering constraint changed the optimum by {a2.objective - a1.objective:.2f}.")
print(" It is an INVALID cut here: it removed the true optimum.")
print()
print("B) Genuinely interchangeable facilities (identical shipping costs)")
b1 = build_facility_model("identical_no_sym", ship_identical, False).solve()
b2 = build_facility_model("identical_sym", ship_identical, True).solve()
print(f" without ordering constraint: obj={b1.objective:.2f} nodes={b1.node_count}")
print(f" with ordering constraint: obj={b2.objective:.2f} nodes={b2.node_count}")
print(f" --> optimum preserved: {abs(b1.objective - b2.objective) < 1e-6}. "
"Here the ordering IS a symmetry break.")
A) Facilities that are NOT interchangeable (shipping costs differ)
without ordering constraint: obj=698.00 nodes=1
with ordering constraint: obj=782.00 nodes=1
--> the ordering constraint changed the optimum by 84.00.
It is an INVALID cut here: it removed the true optimum.
B) Genuinely interchangeable facilities (identical shipping costs)
without ordering constraint: obj=1130.00 nodes=1
with ordering constraint: obj=1130.00 nodes=1
--> optimum preserved: True. Here the ordering IS a symmetry break.
Read case A again. Both solves returned optimal with a zero gap. The second
one is simply wrong about the problem you meant to pose — nothing in the solver
output flags it, because an invalid cut is indistinguishable from a constraint you
genuinely wanted. Symmetry breaking is one of the few modeling tricks that can
silently cost you the optimum, so the interchangeability argument has to be made
before the constraint is written, not after.
In case B the optimum is preserved, which is the validity condition. Note that on an instance this small both variants finish at the root node, so there is no node saving to observe here — presolve already closes it. The payoff from symmetry breaking shows up on instances large enough for B&B to actually enumerate the equivalent permutations, where the search would otherwise revisit \(k!\) relabelings of the same solution.
Modern MILP solvers also detect some symmetries automatically, but the detection is necessarily conservative: it can only exploit symmetries it can prove, so an explicit constraint is often stronger — when you are entitled to it.
7. SOS1 and SOS2 Constraints#
SOS Type 1 (SOS1): At most one variable in the set can be nonzero.
SOS Type 2 (SOS2): At most two consecutive variables can be nonzero (used for piecewise-linear approximations).
SOS constraints are useful when:
Modeling mutually exclusive choices without binary variables
Building piecewise-linear approximations of nonlinear functions
The solver can branch on SOS sets more intelligently than on individual binaries
Example: Mutually exclusive production modes#
A factory can operate in one of three modes (low, medium, high), each with different production and cost levels. Only one mode can be active.
m = dm.Model("sos1_example")
# Production level in each mode (at most one can be nonzero)
x_low = m.continuous("x_low", lb=0, ub=100)
x_med = m.continuous("x_med", lb=0, ub=200)
x_high = m.continuous("x_high", lb=0, ub=500)
# SOS1: at most one mode active
m.sos1([x_low, x_med, x_high], name="mode_selection")
# Revenue minus cost (mode-dependent unit costs: low=5, med=3, high=2)
# Maximize profit
revenue_per_unit = 10.0
m.maximize(
(revenue_per_unit - 5) * x_low
+ (revenue_per_unit - 3) * x_med
+ (revenue_per_unit - 2) * x_high
)
# Resource constraint
m.subject_to(2 * x_low + 3 * x_med + 4 * x_high <= 1000, name="resource")
result = m.solve()
print(f"Status: {result.status}")
print(f"Profit: ${result.objective:.2f}")
print(f"x_low: {result.value(x_low):.2f}")
print(f"x_med: {result.value(x_med):.2f}")
print(f"x_high: {result.value(x_high):.2f}")
print("\nNote: SOS1 ensures at most one variable is nonzero.")
Status: optimal
Profit: $2000.00
x_low: 0.00
x_med: 0.00
x_high: 250.00
Note: SOS1 ensures at most one variable is nonzero.
8. Exercise: Bin Packing#
Given \(n\) items with sizes \(s_i\) and \(B\) bins each with capacity \(C\), find the minimum number of bins needed to pack all items.
Formulation:
subject to: $\(\sum_{b=1}^{B} x_{ib} = 1 \quad \forall i \quad \text{(each item in exactly one bin)}\)\( \)\(\sum_{i=1}^{n} s_i x_{ib} \leq C \cdot u_b \quad \forall b \quad \text{(bin capacity)}\)\( \)\(x_{ib} \in \{0, 1\}, \; u_b \in \{0, 1\}\)$
where \(u_b = 1\) if bin \(b\) is used, and \(x_{ib} = 1\) if item \(i\) is placed in bin \(b\).
Your task: Complete the model below.
# Problem data
item_sizes = np.array([7, 5, 3, 4, 6, 2, 8, 1, 5, 4])
n_items_bp = len(item_sizes)
bin_capacity = 10
n_bins = 6 # upper bound on number of bins
print(f"Items: {n_items_bp}, sizes: {item_sizes}")
print(f"Total size: {item_sizes.sum()}, bin capacity: {bin_capacity}")
print(f"Lower bound on bins: {int(np.ceil(item_sizes.sum() / bin_capacity))}")
# --- YOUR CODE HERE ---
# m = dm.Model("bin_packing")
# u = m.binary("use_bin", shape=(n_bins,)) # bin used?
# x = m.binary("assign", shape=(n_items_bp, n_bins)) # item i in bin b?
#
# m.minimize(dm.sum(u))
#
# # Each item in exactly one bin
# for i in range(n_items_bp):
# m.subject_to(dm.sum(lambda b: x[i, b], over=range(n_bins)) == 1)
#
# # Bin capacity
# for b in range(n_bins):
# m.subject_to(
# dm.sum(lambda i: item_sizes[i] * x[i, b], over=range(n_items_bp)) <= bin_capacity * u[b]
# )
#
# result = m.solve()
# print(f"\nBins needed: {result.objective:.0f}")
Items: 10, sizes: [7 5 3 4 6 2 8 1 5 4]
Total size: 45, bin capacity: 10
Lower bound on bins: 5
# Solution
m = dm.Model("bin_packing")
u = m.binary("use_bin", shape=(n_bins,))
x = m.binary("assign", shape=(n_items_bp, n_bins))
# Minimize number of bins used
m.minimize(dm.sum(u))
# Each item in exactly one bin
for i in range(n_items_bp):
m.subject_to(dm.sum(lambda b: x[i, b], over=range(n_bins)) == 1, name=f"assign_{i}")
# Bin capacity (linking with u)
for b in range(n_bins):
m.subject_to(
dm.sum(lambda i: item_sizes[i] * x[i, b], over=range(n_items_bp)) <= bin_capacity * u[b],
name=f"capacity_{b}",
)
# Symmetry breaking: prefer lower-indexed bins
for b in range(n_bins - 1):
m.subject_to(u[b] >= u[b + 1], name=f"sym_{b}")
result = m.solve(time_limit=60)
print(f"Status: {result.status}")
print(f"Bins needed: {result.objective:.0f}")
print(f"Nodes: {result.node_count}")
print(f"Time: {result.wall_time:.3f}s")
print()
# Show assignment
assign = result.value(x)
for b in range(n_bins):
items_in_bin = [i for i in range(n_items_bp) if assign[i, b] > 0.5]
if items_in_bin:
sizes = [item_sizes[i] for i in items_in_bin]
print(f"Bin {b}: items {items_in_bin}, sizes {sizes}, total={sum(sizes)}")
Status: optimal
Bins needed: 5
Nodes: 1
Time: 0.013s
Bin 0: items [1, 8], sizes [np.int64(5), np.int64(5)], total=10
Bin 1: items [0], sizes [np.int64(7)], total=7
Bin 2: items [6], sizes [np.int64(8)], total=8
Bin 3: items [2, 3, 5, 7], sizes [np.int64(3), np.int64(4), np.int64(2), np.int64(1)], total=10
Bin 4: items [4, 9], sizes [np.int64(6), np.int64(4)], total=10
9. Summary#
In this tutorial we covered:
Topic |
Key takeaway |
|---|---|
LP to MILP |
Integrality constraints make LPs NP-hard; Branch-and-Bound is the workhorse algorithm |
B&B mechanics |
LP relaxation gives lower bounds; incumbent gives upper bound; prune when LB \(\geq\) UB |
Knapsack |
Classic 0-1 MILP; LP relaxation is typically tight |
Facility location |
Big-M vs indicator constraints; formulation strength matters |
Job-shop scheduling |
Disjunctive constraints for sequencing; big-M for “either/or” logic |
Symmetry breaking |
Ordering constraints eliminate equivalent solutions and reduce nodes |
SOS1 constraints |
Model mutually exclusive choices without introducing binary variables |
Bin packing |
Combining assignment and capacity constraints |
Practical tips for formulating MILPs#
Start with the LP relaxation. If it is weak, tighten it with valid inequalities or reformulation.
Use the smallest possible big-M. Overly large \(M\) values weaken the LP relaxation.
Break symmetries whenever you have interchangeable objects.
Check the gap. If gap = 0, the solution is provably optimal. If not, consider running longer or reformulating.
Node count is your friend. It tells you how hard the instance was for B&B.
Next steps#
Mixed-Integer Quadratic Programming (MIQP): Adding quadratic objectives to MILP, covered in the MIQP tutorial notebook.
MINLP: Full nonlinear constraints, covered in
notebooks/minlp_examples.ipynb.Advanced B&B features: Cutting planes, learned relaxations — see
notebooks/advanced_features.ipynb.
References#
[Land and Doig, 1960] — Branch and Bound algorithm
[Dakin, 1965] — B&B refinements
[Wolsey, 1998] — Integer Programming (textbook)
[Nemhauser and Wolsey, 1988] — Integer and Combinatorial Optimization (textbook)
[Bertsimas and Tsitsiklis, 1997] — Introduction to Linear Optimization