Cutting Stock

Cutting Stock

You are an expert in the one-dimensional cutting stock problem (1D-CSP): the canonical pattern-based optimization problem and the original application of column generation (Gilmore & Gomory, 1961, "A linear programming approach to the cutting-stock problem"). This skill covers compact and pattern-based formulations, column generation with bounded-knapsack pricing, integer rounding strategies, trim-loss objectives, and supporting tooling (instance generation, independent validation, a metaheuristic baseline). Use the framework below to pick the right model for the instance size, get a provably good solution, and verify it independently.

Initial Assessment

Establish these facts before proposing a model or writing code:

• Instance dimensions. Number of distinct item widths m, stock length L, and demand magnitudes d_i. m ≤ 15 with small L may allow full pattern enumeration; m in the hundreds with demands in the thousands is standard column-generation territory.
• Width data type. Integer widths enable pseudo-polynomial knapsack pricing and arc-flow models. Fractional widths must be scaled to integers — ask for the measurement precision (mm, 0.1 mm) and check the scaled L stays manageable.
• Stock assortment. One stock length or several? Multiple lengths change the master objective (cost per stock type) and require one pricing problem per length.
• Objective. Minimize number of rolls, total trim loss, or material cost? With identical rolls these align, but only under a stated overproduction policy — confirm whether cutting more pieces than demanded is waste, usable inventory, or forbidden.
• Demand semantics. Cover-at-least (≥ d_i, the default, gives nonnegative duals and a clean pricing problem) or meet-exactly (= d_i, harder: duals can be negative, master can be infeasible with few columns)?
• Side constraints. Maximum number of distinct patterns (setup costs), maximum pieces per pattern (knife count), pattern run-length limits, due dates per order. These decide between vanilla Gilmore-Gomory and an extension.
• Optimality requirement. Is ceil(LP bound) or +1 roll acceptable (almost always reached by rounding heuristics), or is a proven optimum required (branch-and-price or arc-flow)?
• Solver availability. Gurobi licensed? If not, the same pattern applies with HiGHS/CBC for the master and a hand-written DP for pricing — only the master LP/IP calls change.
• Time budget. Column generation on m ≤ 200 converges in seconds; a compact MIP on the same instance may not finish in hours because of symmetry.
• Validation path. Agree up front on an independent feasibility checker (pattern widths, demand coverage) so the model and the check do not share code.

Problem Variants and Formulation

Formal definition

Given stock pieces of integer length $L$ and item types $i \in I = \{1,\dots,m\}$ with width $wi \in \mathbb{Z}{>0}$, $wi \le L$, and demand $di \in \mathbb{Z}{>0}$: cut all demanded items from the minimum number of stock pieces. A cutting pattern is a vector $a = (a1,\dots,am) \in \mathbb{Z}{\ge 0}^m$ with $\sumi wi ai \le L$; $ai$ counts copies of item $i$ cut from one stock piece.

Pattern-based (Gilmore-Gomory) model. Let $J$ index all feasible patterns and $x_j$ the number of stock pieces cut with pattern $j$:

$$\min \sum{j \in J} xj \quad \text{s.t.} \quad \sum{j \in J} a{ij}\, xj \ \ge\ di \quad \forall i \in I, \qquad xj \in \mathbb{Z}{\ge 0}.$$

$|J|$ grows exponentially in $m$, so the LP relaxation is solved by column generation: the pricing problem $\max\{\sumi \pii ai : \sumi wi ai \le L,\ a \in \mathbb{Z}{\ge 0}^m\}$ is a bounded integer knapsack over the dual prices $\pii$ of the demand rows. A pattern prices out (improves the master) iff its reduced cost $1 - \sumi \pii ai list[list[int]]: """FFD on the expanded item list. Returns rolls as lists of item-type indices.""" items = [i for i, d in enumerate(demands) for in range(d)] items.sort(key=lambda i: -widths[i]) rolls: list[list[int]] = [] space: list[int] = [] for i in items: for r in range(len(rolls)): if widths[i] None: """Each item type i is cut at least d_i times across all rolls.""" for i in range(len(data["widths"])): model.addConstr( gp.quicksum(z[i, k] for k in range(data["K"])) >= data["demands"][i], name=f"demand[{i}]", )

def addcapacityconstraints(model: gp.Model, y: gp.tupledict, z: gp.tupledict, data: dict) -> None: """Width cut from roll k fits in the stock length; links z to the roll-open var y.""" m = len(data["widths"]) for k in range(data["K"]): model.addConstr( gp.quicksum(data["widths"][i] * z[i, k] for i in range(m)) None: """Order the rolls: roll k+1 may be used only if roll k is used.""" for k in range(data["K"] - 1): model.addConstr(y[k] >= y[k + 1], name=f"sym[{k}]")

def solvekantorovich(widths: list[int], demands: list[int], capacity: int, timelimit: float = 60.0) -> tuple[int, list[list[int]]]: """Compact MIP for 1D cutting stock. Returns (rollsused, per-roll patterns).""" K = len(firstfitdecreasing(widths, demands, capacity)) # valid upper bound data = {"widths": widths, "demands": demands, "capacity": capacity, "K": K} m = len(widths) model = gp.Model("cspkantorovich") model.Params.OutputFlag = 0 model.Params.TimeLimit = timelimit y = model.addVars(K, vtype=GRB.BINARY, name="y") z = model.addVars(m, K, vtype=GRB.INTEGER, lb=0, name="z") adddemandconstraints(model, z, data) addcapacityconstraints(model, y, z, data) addsymmetrybreakingconstraints(model, y, data) 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 solution, status {model.Status}") used = [k for k in range(K) if y[k].X > 0.5] patterns = [[round(z[i, k].X) for i in range(m)] for k in used] return len(used), patterns

widths, demands, capacity = [6, 5, 4, 3], [2, 2, 3, 4], 10 rolls, patterns = solve_kantorovich(widths, demands, capacity) print(rolls, patterns)

Expected: rolls = 5. Material bound ceil(46/10) = 5, so 5 is provably optimal.

Even with the `y[k] >= y[k+1]` ordering, the compact model degrades quickly: identical rolls remain interchangeable in the `z` block, and the LP bound is the material bound. Treat it as a correctness reference for `m ≲ 10` and total demand ≲ 50, not as a production model.

### Pattern enumeration + pattern IP (small instances)

A pattern is **maximal** if no further item fits into its residual width. Restricting to maximal patterns never hurts the roll-count objective with `≥` demand rows.

```python
import gurobipy as gp
from gurobipy import GRB

def enumerate_maximal_patterns(widths: list[int], capacity: int) -> list[list[int]]:
    """DFS over all maximal cutting patterns (no item fits in the residual width)."""
    m = len(widths)
    patterns: list[list[int]] = []

    def extend(i: int, residual: int, current: list[int]) -> None:
        if i == m:
            if all(residual  tuple[int, dict[int, int]]:
    """Set-covering style pattern IP over an explicit pattern list."""
    model = gp.Model("csp_pattern_ip")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    x = model.addVars(len(patterns), vtype=GRB.INTEGER, lb=0, name="x")
    for i in range(len(demands)):
        model.addConstr(
            gp.quicksum(patterns[j][i] * x[j] for j in range(len(patterns)))
            >= demands[i],
            name=f"cover[{i}]",
        )
    model.setObjective(x.sum(), GRB.MINIMIZE)
    model.optimize()
    if model.Status not in (GRB.OPTIMAL, GRB.TIME_LIMIT) or model.SolCount == 0:
        raise RuntimeError(f"no solution, status {model.Status}")
    use = {j: round(x[j].X) for j in range(len(patterns)) if x[j].X > 0.5}
    return round(model.ObjVal), use

widths, demands, capacity = [6, 5, 4, 3], [2, 2, 3, 4], 10
pats = enumerate_maximal_patterns(widths, capacity)
obj, use = solve_pattern_ip(pats, demands)
print(len(pats), obj, {tuple(pats[j]): c for j, c in use.items()})
# Expected: 8 maximal patterns; optimum 5 rolls,
# e.g. (1,0,1,0) x2, (0,2,0,0) x1, (0,0,1,2) x2.

When enumeration is feasible this is the simplest provably exact method: the IP sees every column, so there is no integrality question beyond the IP solve itself.

Gilmore-Gomory Column Generation

The workhorse. The restricted master LP holds a subset of patterns; bounded-knapsack pricing over the duals generates improving patterns until none has negative reduced cost. The loop below adds columns incrementally with gp.Column instead of rebuilding the master. For stabilization, convergence theory, and branch-and-price mechanics see the column-generation skill; for the pricing DP family (bounded/unbounded knapsack, branch-and-bound alternatives) see the knapsack-problems skill.

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

def solve_pricing(widths: np.ndarray, duals: np.ndarray, capacity: int,
                  max_copies: np.ndarray) -> tuple[np.ndarray, float]:
    """Bounded integer knapsack by DP: max duals @ a, s.t. widths @ a  new[c] + 1e-12:
                    new[c] = v
                    take[i, c] = k
        dp = new
    pattern = np.zeros(m, dtype=np.int64)
    c = capacity
    for i in range(m - 1, -1, -1):
        pattern[i] = take[i, c]
        c -= int(take[i, c] * widths[i])
    return pattern, float(dp[capacity])

def column_generation(widths: list[int], demands: list[int], capacity: int,
                      tol: float = 1e-9, max_iter: int = 1000
                      ) -> tuple[list[list[int]], float, int]:
    """Gilmore-Gomory CG for the master LP. Returns (patterns, z_LP, iterations)."""
    m = len(widths)
    patterns = []
    for i in range(m):                       # homogeneous starting columns
        col = [0] * m
        col[i] = capacity // widths[i]
        patterns.append(col)
    model = gp.Model("csp_master_lp")
    model.Params.OutputFlag = 0
    cover = [model.addConstr(gp.LinExpr() >= demands[i], name=f"cover[{i}]")
             for i in range(m)]
    for j, col in enumerate(patterns):
        nz = [i for i in range(m) if col[i]]
        model.addVar(obj=1.0, lb=0.0, name=f"x[{j}]",
                     column=gp.Column([float(col[i]) for i in nz],
                                      [cover[i] for i in nz]))
    model.ModelSense = GRB.MINIMIZE
    w = np.asarray(widths, dtype=np.int64)
    d = np.asarray(demands, dtype=np.int64)
    for it in range(1, max_iter + 1):
        model.optimize()
        if model.Status != GRB.OPTIMAL:
            raise RuntimeError(f"master LP status {model.Status}")
        duals = np.array([c.Pi for c in cover])
        pattern, value = solve_pricing(w, duals, capacity, d)
        if value = -tol
            return patterns, model.ObjVal, it
        patterns.append(pattern.tolist())
        nz = [i for i in range(m) if pattern[i]]
        model.addVar(obj=1.0, lb=0.0, name=f"x[{len(patterns) - 1}]",
                     column=gp.Column([float(pattern[i]) for i in nz],
                                      [cover[i] for i in nz]))
    raise RuntimeError("column generation hit max_iter without converging")

def solve_ip_over_columns(patterns: list[list[int]], demands: list[int],
                          time_limit: float = 30.0) -> tuple[int, list[int]]:
    """Restricted-master IP: integer solve over the generated columns only."""
    model = gp.Model("csp_master_ip")
    model.Params.OutputFlag = 0
    model.Params.TimeLimit = time_limit
    x = model.addVars(len(patterns), vtype=GRB.INTEGER, lb=0, name="x")
    for i in range(len(demands)):
        model.addConstr(
            gp.quicksum(patterns[j][i] * x[j] for j in range(len(patterns)))
            >= demands[i],
            name=f"cover[{i}]",
        )
    model.setObjective(x.sum(), GRB.MINIMIZE)
    model.optimize()
    if model.Status not in (GRB.OPTIMAL, GRB.TIME_LIMIT) or model.SolCount == 0:
        raise RuntimeError(f"no solution, status {model.Status}")
    return round(model.ObjVal), [round(x[j].X) for j in range(len(patterns))]

widths, demands, capacity = [6, 5, 4, 3], [2, 2, 3, 4], 10
patterns, z_lp, iters = column_generation(widths, demands, capacity)
z_ip, counts = solve_ip_over_columns(patterns, demands)
print(f"z_LP={z_lp:.4f} ceil={int(np.ceil(z_lp - 1e-9))} z_IP={z_ip} iters={iters}")
# Expected: z_LP about 4.67, ceil(z_LP) = 5, z_IP = 5 -> proven optimal
# (material bound is 4.6; the LP bound dominates it).

Three facts to internalize:

1. The stopping test is on the pricing optimum, not on master progress. The master objective can stall for many iterations (degeneracy) while the duals still move; stop only when $\maxa \sumi \pii ai tuple[list[list[int]], list[int]]:

"""Round the LP master solution down; cover residual demand by FFD. Returns (allpatterns, counts) including the repair patterns.""" a = np.asarray(patterns, dtype=np.int64) xdn = np.floor(np.asarray(xlp) + 1e-9).astype(np.int64) residual = np.maximum(np.asarray(demands) - a.T @ xdn, 0) items = np.repeat(np.arange(len(widths)), residual) # expand residual items order = sorted(range(len(items)), key=lambda r: -widths[items[r]]) rolls: list[list[int]] = [] loads: list[int] = [] for r in order: t = int(items[r]) for k in range(len(rolls)): if loads[k] + widths[t] optimal.

Report whichever of {IP over columns, round-down + repair} is better, together with $\lceil z_{LP} \rceil$; if they meet the bound, say "optimal", otherwise report the absolute gap in rolls (it is almost always 0 or 1 — see MIRUP above).

## Metaheuristic Baseline: Random-Key GA

A metaheuristic is the right tool when side constraints break the knapsack pricing structure (knife limits, sequencing, pattern-dependent costs) or when no LP solver is available. For plain 1D-CSP it is a baseline, not a competitor to CG. A biased random-key GA with a first-fit decoder fits well: the decoder maps any key vector to a feasible solution, so no constraint handling is needed. Framework details, elite/mutant proportions, and decoder design guidance are in the **biased-random-key-genetic-algorithm** skill; the grouping-fitness idea is from Falkenauer (1996), "A hybrid grouping genetic algorithm for bin packing".

```python
import numpy as np
from collections import Counter

def decode_first_fit(keys: np.ndarray, item_types: np.ndarray, widths: list[int],
                     capacity: int) -> tuple[list[list[int]], list[int]]:
    """Sort items by key, first-fit into rolls. Returns (rolls as type lists, loads)."""
    rolls: list[list[int]] = []
    loads: list[int] = []
    for idx in np.argsort(keys):
        t = int(item_types[idx])
        for r in range(len(rolls)):
            if loads[r] + widths[t]  float:
    """Falkenauer-style: roll count minus mean squared fill (prefer full rolls)."""
    fills = np.asarray(loads, dtype=float) / capacity
    return len(loads) - float(np.mean(fills**2))

def brkga_cutting_stock(widths: list[int], demands: list[int], capacity: int,
                        pop_size: int = 60, n_gen: int = 200, elite: float = 0.2,
                        mutant: float = 0.15, rho: float = 0.7, seed: int = 0
                        ) -> tuple[int, Counter]:
    """BRKGA over item-permutation random keys with a first-fit decoder."""
    rng = np.random.default_rng(seed)
    item_types = np.repeat(np.arange(len(widths)), demands)
    n = len(item_types)
    n_e = max(1, int(elite * pop_size))
    n_m = max(1, int(mutant * pop_size))
    n_o = pop_s

…

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