Bin Packing

Bin Packing

You are an expert in one-dimensional bin packing and its standard variants. This skill covers construction heuristics with worst-case guarantees (FF, BF, FFD, BFD), the L1 and L2 lower bounds, the compact assignment MIP with symmetry breaking, the arc-flow exact model, item conflicts, variable bin sizes, and the relation to cutting stock. Use the framework below to choose the right bound, heuristic, and exact model for the instance at hand, and to validate every solution independently.

Initial Assessment

Establish the following before proposing a method.

• Instance size. Get n (number of items), the capacity C, and the number of distinct

item sizes d. If d C makes the instance infeasible; items with w_i = 0 should be removed before modeling.

• Variant detection. Identical bins, or multiple bin types with different capacities and

costs (variable-sized bin packing)? Are there incompatible item pairs that may not share a bin (bin packing with conflicts)? Cardinality limits per bin?

• Objective check. Confirm the goal is minimizing the number of bins. If the number of bins

is fixed and the goal is balancing loads, that is multiprocessor scheduling (P||Cmax), a different problem with different methods.

• Online vs offline. Do all items arrive up front? Online arrival changes the achievable

guarantees (no algorithm beats competitive ratio ~1.54) and rules out sorting-based methods.

• Exactness need. Is a provably optimal count required, or is "within one bin of a lower

bound, certified" enough? FFD plus L2 often closes the gap without any solver.

• Time budget and solver access. Seconds or hours? Is a Gurobi license available, or should

the model run on HiGHS / CP-SAT? Arc-flow models can be large: estimate C * d arcs first.

• Solution artifact. Does the user need only the bin count, or the explicit item-to-bin

assignment? The assignment matters for downstream use and for validation.

• Scale of repetition. One instance, or thousands solved in a loop (e.g., inside a pricing

or simulation routine)? Repetition pushes toward O(n log n) heuristics with cached bounds.

Problem Definition, Bounds, and Model Choice

Formal definition. Given items i in I = {1, ..., n} with weights w_i > 0 and identical bins of capacity C >= max_i w_i, partition I into the minimum number of subsets (bins) such that the total weight in each subset is at most C. The compact (Kantorovich) formulation over a bin index set B = {1, ..., U} (U = any valid upper bound):

$$ \min \sum{b \in B} yb \quad \text{s.t.} \quad \sum{b \in B} x{ib} = 1 \;\; \forall i \in I, \qquad \sum{i \in I} wi\, x{ib} \le C\, yb \;\; \forall b \in B, \qquad x{ib},\, yb \in \{0, 1\}. $$

The pattern (set-covering / Gilmore-Gomory) formulation uses the set P of all maximal feasible packings of a single bin:

$$ \min \sum{p \in P} \lambdap \quad \text{s.t.} \quad \sum{p \ni i} \lambdap \ge 1 \;\; \forall i \in I, \qquad \lambdap \in \mathbb{Z}{\ge 0}. $$

Its LP relaxation z_LP is far stronger than the compact LP (which only gives the trivial L1 bound) and is solved by column generation with a knapsack pricing problem.

Complexity. Bin packing is strongly NP-hard (reduction from 3-Partition; Garey & Johnson 1979, "Computers and Intractability"). No polynomial algorithm approximates within a factor below 3/2 unless P = NP (deciding 2 vs 3 bins encodes PARTITION). There is an APTAS (Fernandez de la Vega & Lueker 1981) and an additive OPT + O(log OPT) algorithm via the pattern LP (Hoberg & Rothvoss 2017). In practice, instances with thousands of items are routinely solved to optimality by arc-flow or branch-and-price.

Lower bounds.

• L1 = ceil(sum_i w_i / C): the continuous bound, equal to the compact LP relaxation. Can be

as weak as OPT/2 (e.g., all items slightly larger than C/2).

• L2 (Martello & Toth 1990, "Knapsack Problems"): for a threshold alpha in [0, C/2], count

items larger than C - alpha (one bin each) plus items in (C/2, C - alpha], then charge the items in [alpha, C/2] against the residual space of the latter group. Maximize over alpha. Computable in O(n log n); worst-case ratio 2/3 of OPT, and frequently tight.

• ceil(z_LP) from the pattern LP: the strongest practical bound. The MIRUP conjecture

(Scheithauer & Terno 1995) states OPT list[list[int]]: """Pack items in order (default: input order), each into the first bin that fits.""" idx = list(order) if order is not None else list(range(len(weights))) bins: list[list[int]] = [] loads: list[float] = [] for i in idx: for b, load in enumerate(loads): if load + weights[i] list[list[int]]: """Pack items in order`, each into the feasible bin with least residual space.""" idx = list(order) if order is not None else list(range(len(weights))) bins: list[list[int]] = [] loads: list[float] = [] for i in idx: bestb, bestresid = -1, float("inf") for b, load in enumerate(loads): resid = capacity - load - weights[i] if 0 = 0: bins[bestb].append(i) loads[bestb] += weights[i] else: bins.append([i]) loads.append(weights[i]) return bins

def decreasing_order(weights: Sequence[float]) -> list[int]: """Indices sorted by non-increasing weight (the 'D' in FFD/BFD).""" return sorted(range(len(weights)), key=lambda i: -weights[i])

weights = [7, 6, 5, 4, 3, 2, 1] ffdbins = firstfit(weights, capacity=10, order=decreasingorder(weights)) bfdbins = bestfit(weights, capacity=10, order=decreasingorder(weights)) print(len(ffdbins), ffdbins)

Expected: 3 [[0, 4], [1, 3], [2, 5, 6]] — loads 10, 10, 8; L1 = ceil(28/10) = 3,

so FFD is provably optimal here without any solver.

Always compute both lower bounds before judging a heuristic solution. If the heuristic count
equals `max(L1, L2)`, you are done.

```python
import math

import numpy as np

def lower_bound_l1(weights: np.ndarray, capacity: int) -> int:
    """Continuous bound: ceil(total weight / capacity). Equals the compact LP bound."""
    return math.ceil(float(weights.sum()) / capacity)

def lower_bound_l2(weights: np.ndarray, capacity: int) -> int:
    """Martello-Toth L2 bound, O(n log n) over candidate thresholds.

    For threshold alpha: items > C - alpha each need their own bin (J1); items in
    (C/2, C - alpha] need their own bin too (J2) but keep usable residual space;
    items in [alpha, C/2] (J3) are charged against that residual space. It suffices
    to scan alpha over 0 and the distinct weights  capacity - alpha]
        j2 = w[(w > capacity / 2) & (w = alpha) & (w = OPT/2` and `L2 >= (2/3) OPT`
asymptotically (Martello & Toth 1990). When `L2  list[tuple[int, int]]:
    """Canonical symmetry reduction: item i may only use bins 0..min(i, n_bins-1)."""
    return [(i, b) for i in range(n_items) for b in range(min(i + 1, n_bins))]

def add_assignment_constraints(model: gp.Model, x: gp.tupledict,
                               n_items: int, n_bins: int) -> None:
    """Each item is placed in exactly one of its allowed bins."""
    for i in range(n_items):
        model.addConstr(
            gp.quicksum(x[i, b] for b in range(min(i + 1, n_bins))) == 1,
            name=f"assign[{i}]")

def add_capacity_constraints(model: gp.Model, x: gp.tupledict, y: gp.tupledict,
                             weights: np.ndarray, capacity: float,
                             n_bins: int) -> None:
    """Bin load at most capacity, and only if the bin is open (links x to y)."""
    n_items = len(weights)
    for b in range(n_bins):
        model.addConstr(
            gp.quicksum(float(weights[i]) * x[i, b] for i in range(b, n_items))
             None:
    """Open bins form a prefix: bin b+1 may open only if bin b is open."""
    for b in range(n_bins - 1):
        model.addConstr(y[b] >= y[b + 1], name=f"prefix[{b}]")

def add_conflict_constraints(model: gp.Model, x: gp.tupledict,
                             conflicts: list[tuple[int, int]],
                             n_bins: int) -> None:
    """Conflicting items may not share a bin (only bins where both can appear)."""
    for i, j in conflicts:
        i, j = min(i, j), max(i, j)
        for b in range(min(i + 1, n_bins)):
            model.addConstr(x[i, b] + x[j, b]  tuple[list[list[int]], float]:
    """Compact assignment MIP. Pass n_bins = FFD count as the upper bound."""
    n_items = len(weights)
    if n_bins is None:
        n_bins = n_items
    model = gp.Model("bin_packing")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    x = model.addVars(allowed_pairs(n_items, n_bins), vtype=GRB.BINARY, name="x")
    y = model.addVars(n_bins, vtype=GRB.BINARY, name="y")
    add_assignment_constraints(model, x, n_items, n_bins)
    add_capacity_constraints(model, x, y, np.asarray(weights), capacity, n_bins)
    add_symmetry_breaking_constraints(model, y, n_bins)
    if conflicts:
        add_conflict_constraints(model, x, conflicts, n_bins)
    model.setObjective(y.sum(), GRB.MINIMIZE)
    model.optimize()
    if model.Status not in (GRB.OPTIMAL, GRB.TIME_LIMIT) or model.SolCount == 0:
        raise RuntimeError(f"no feasible solution found, status {model.Status}")
    bins = [[i for i in range(n_items) if (i, b) in x and x[i, b].X > 0.5]
            for b in range(n_bins)]
    return [items for items in bins if items], model.ObjBound

w = np.array([7, 6, 5, 4, 3, 2, 1])
bins, bound = solve_bin_packing_mip(w, capacity=10, n_bins=4)
print(len(bins), round(bound))
# Expected: 3 3 — three bins, and ObjBound = 3 proves optimality.

Expect the compact model to handle up to roughly 100-200 items with tight bounds; beyond that, symmetry and the weak LP make it stall. Switch to arc-flow or column generation. For general gurobipy modeling idioms (parameters, status handling, warm starts), see milp-modeling-gurobi.

Arc-flow model

The arc-flow formulation (Valério de Carvalho 1999, "Exact solution of bin-packing problems using column generation and branch-and-bound") models one bin as a path in a graph whose nodes are partial loads 0..C. An item arc (u, u + w) packs one item of size w on top of load u; loss arcs (u, u + 1) absorb unused capacity. The total flow value equals the number of bins, and demand constraints force each size class to be used the right number of times. Its LP relaxation is equal in strength to the Gilmore-Gomory pattern LP, so the model usually solves at the root node. Size is pseudo-polynomial — O(C·d) arcs — so it needs integer weights and a moderate capacity. The version below restricts arc tails to subset-sum-reachable nodes, the basic graph reduction; see Advanced Techniques for stronger reductions.

from collections import defaultdict

import gurobipy as gp
import numpy as np
from gurobipy import GRB

def reachable_nodes(sizes: list[int], counts: list[int], capacity: int) -> list[int]:
    """Subset-sum DP respecting multiplicities: only these loads can be arc tails."""
    reach = np.zeros(capacity + 1, dtype=bool)
    reach[0] = True
    for w, c in zip(sizes, counts):
        for _ in range(c):
            new = reach.copy()
            new[w:] |= reach[:-w]
            if (new == reach).all():
                break
            reach = new
    return [u for u in range(capacity + 1) if reach[u]]

def solve_arc_flow(weights: np.ndarray, capacity: int,
                   time_limit: float = 60.0) -> int:
    """Arc-flow bin packing model; returns the optimal number of bins."""
    sizes_arr, counts_arr = np.unique(np.asarray(weights, dtype=int),
                                      return_counts=True)
    sizes = [int(s) for s in sizes_arr]
    counts = [int(c) for c in counts_arr]
    nodes = reachable_nodes(sizes, counts, capacity)
    item_arcs = [(u, u + w, w) for w in sizes for u in nodes if u + w = 1 else 0)
        outflow = gp.quicksum(f[a] for a in arcs_out[u]) \
            + (loss[u] if u  ceil = 3, so the root node closes the instance.

Recovering the explicit packing from an arc-flow solution means decomposing the flow into z source-to-sink paths (each path is one bin); do this greedily on the integer flow values. For demand-aggregated data and trim-loss objectives, the same machinery is developed further in cutting-stock; the pricing view of the same LP is in column-generation.

Metaheuristic: Random-Key GA with First-Fit Decoder

When FFD leaves a gap of 2+ bins and exact models are too slow (huge n, fractional weights, messy side constraints), a biased random-key GA over packing orders is a robust choice: keys in [0,1)^n are decoded by sorting into an item order and packing with first fit. Every decoded solution is feasible by construction, and the only problem-specific code is the decoder. The fitness pairs the bin count with Falkenauer's fill measure mean((load_b / C)^2) (Falkenauer 1996, "A hybrid grouping genetic algorithm for bin packing") so that, among equal bin counts, packings with fuller bins — which are closer to dropping a bin — rank higher. Population operations are vectorized; see the biased-random-key-genetic-algorithm skill for the framework details (elite/mutant partitioning, biased crossover) and hyper-heuristics for managing a pool of packing rules instead of a single decoder.

import numpy as np

def first_fit_decode(perm: np.ndarray, weights: np.ndarray,
                     capacity: float) -> tuple[int, float, list[list[int]]]:
    """Decode an item permutation with first fit.

    Returns (bin count, fill measure mean((load/C)^2), bins). The fill measure
    breaks ties between equal bin counts in favor of fuller bins.
    """
    bins: list[list[int]] = []
    loads: list[float] = []
    for i in perm:
        w = float(weights[i])
        for b, load in enumerate(loads):
            if load + w  tuple[int, list[list[int]]]:
    """BRKGA-style random-key GA for bin packing; returns (bin count, bins)."""
    rng = np.random.default_rng(seed)
    n = len(weights)
    n_elite = max(1, int(elite_frac * pop_size))
    n_mut = max(1, int(mutant_frac * pop_size))
    n_off = pop_size - n_elite - n_mut
    pop = rng.random((pop_size, n))
    pop[0] = np.argsort(np.argsort(-np.asarray(weights, dtype=float))) / n  # FFD seed
    best_count, best_bins = n + 1, [[i] for i in range(n)]
    for _ in range(generations):
        perms = np.argsort(pop, axis=1)
        results = [first_fit_decode(p, weights, capacity) for p in perms]
        scores = np.array([(r[0], r[1]) for r in results])
        rank = np.lexsort((-scores[:, 1], scores[:, 0]))  # bins asc, fill desc
        if int(scores[rank[0], 0])  np.ndarray:
    """Integer weights uniform in [lo*C, hi*C]. Classic hard band: lo=0.2, hi=0.8."""
    rng = np.random.default_rng(seed)
    low = max(1, int(round(lo * capacity)))
    high = max(low, int(round(hi * capacity)))
    return rng.integers(low, high + 1, size=n)

def triplet_instance(n_bins: int, capacity: int, seed: int) -> tuple[np.ndarray, int]:
    """Falkenauer-style triplet instance with known optimum = n_bins.

    Each generated bin holds items (a, b, c) with a + b + c = C exactly,
    a in [0.38C, 0.49C], b and c in [0.25C, 0.5C). Zero slack overall.
    """
    rng = np.random.default_rng(seed)
    a = rng.integers(int(0.38 * capacity), int(0.49 * capacity) + 1, size=n_bins)
    b = rng.integers(int(0.25 * capacity), (capacity - a) // 2 + 1)
    c = capacity - a - b
    weights = np.concatenate([a, b, c])
    rng.shuffle(weights)
    return weights, n_bins

def conflict_graph(n: int, density: float, seed: int) -> list[tuple[int, int]]:
    """Random conflict pairs (i  tuple[bool, list[str]]:
    """Independent feasibility check; shares no code with any solution method.

    Verifies: every item packed exactly once, every bi

…

## Source & license

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

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

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