hajibabaie Dantzig-Wolfe Decomposition
You are an expert in Dantzig-Wolfe (DW) decomposition for linear and mixed-integer programs. This skill covers detecting block-angular structure in a constraint matrix, reformulating the original (compact) model into a master problem with convexity constraints plus independent pricing subproblems, solving the reformulation by column generation, and interpreting the resulting bound against the LP relaxation, the Lagrangian dual, and the integer optimum. Use the framework below to decide whether a model is worth decomposing, to execute the reformulation correctly, and to verify the bound relationships on the user's instance.
Initial Assessment
Establish the following before reformulating anything:
- Structure. Which constraint rows couple otherwise-independent variable groups? Ask the user to name the natural blocks (plants, vehicles, machines, periods, scenarios). If they cannot, run structure detection on the constraint matrix (section below).
- Linking fraction. Count linking rows m0 versus total rows. DW pays off when m0 is small relative to the block rows — a useful rule of thumb is linking rows below 10-20% of all rows.
- Problem class. LP or MIP? If MIP, locate the integrality: integer variables inside blocks make the DW bound potentially stronger than the LP bound; integrality that lives only in the linking rows gains nothing from convexification.
- Block inventory. Number of blocks K, variables per block, and whether the blocks are identical (same costs, same constraint data). Identical blocks call for the aggregated master, which removes symmetry.
- Pricing tractability. What does one block look like in isolation? A knapsack, a shortest path, a small assignment, a small MIP? The whole method stands or falls on solving the pricing problem quickly and repeatedly.
- Boundedness. Are the block polyhedra bounded? If not, the implementation must handle extreme rays, not just extreme points.
- Goal. A tighter dual bound, a faster LP solve on a huge structured model, or an integer optimum? The first two end with column generation; the third requires branch-and-price (hand off to the column-generation skill).
- Baseline. Solve the compact model (or its LP relaxation) first. Record z_LP, the MIP gap, and the time. Without this baseline you cannot say whether DW helped.
- Solver access. Gurobi license for master and pricing? If not, plan for HiGHS as the master LP solver and GCG/SCIP for an automatic end-to-end alternative.
- Stopping policy. Run column generation to proven optimality, or stop early on a Lagrangian-bound gap? Agree on the tolerance up front.
- Time budget. Per-iteration cost is one RMP LP plus K pricing solves. Estimate iterations in the tens-to-hundreds range and check the budget supports that.
Decomposition Anatomy
Block-angular form
DW applies to problems whose constraint matrix is block-angular: a few linking rows across all variables, then independent diagonal blocks.
$$ z_{\mathrm{IP}} ;=; \min ;\sum_{k=1}^{K} c_k^{\top} x_k \quad \text{s.t.} \quad \sum_{k=1}^{K} A_k x_k ;{\le,=,\ge}; b ;;[\pi], \qquad x_k \in X_k, ;; k = 1,\dots,K, $$
where $X_k = {x \in \mathbb{Z}+^{n_k} \text{ (or } \mathbb{R}+^{n_k}\text{)} : D_k x \le d_k}$ collects the block-k constraints and the linking constraints have $m_0$ rows with duals $\pi$.
Reformulation by convexification
By the Minkowski-Weyl theorem, every point of $\operatorname{conv}(X_k)$ is a convex combination of its extreme points ${x_k^p}{p \in P_k}$ plus a conic combination of its extreme rays ${r_k^q}{q \in Q_k}$. Substituting
$$ x_k = \sum_{p \in P_k} \lambda_{kp}, x_k^p + \sum_{q \in Q_k} \mu_{kq}, r_k^q, \qquad \sum_{p \in P_k} \lambda_{kp} = 1, \quad \lambda, \mu \ge 0 $$
into the linking constraints yields the DW master problem:
$$ z_{\mathrm{DW}} = \min \sum_{k,p} (c_k^{\top} x_k^p), \lambda_{kp} + \sum_{k,q} (c_k^{\top} r_k^q), \mu_{kq} $$
$$ \text{s.t.} \quad \sum_{k,p} (A_k x_k^p), \lambda_{kp} + \sum_{k,q} (A_k r_k^q), \mu_{kq} ;{\le,=,\ge}; b ;;[\pi], \qquad \sum_{p} \lambda_{kp} = 1 ;;[\sigma_k], ;; k=1,\dots,K. $$
The per-block equations $\sum_p \lambda_{kp} = 1$ are the convexity constraints; their duals $\sigma_k$ price block membership. The master has one column per extreme point or ray — exponentially many — so it is solved by column generation over a restricted master problem (RMP): solve the RMP, read $(\pi, \sigma)$, and for each block solve the pricing problem
$$ \bar{c}k ;=; \min{x \in X_k} ;(c_k - A_k^{\top}\pi)^{\top} x ;-; \sigma_k . $$
If $\bar{c}_k < 0$, add the optimal point as a column (coefficients $A_k x$, convexity coefficient 1). If pricing is unbounded, the returned extreme ray $r$ with $(c_k - A_k^{\top}\pi)^{\top} r < 0$ enters as a ray column (no convexity coefficient). When no block prices out, the master is optimal (Dantzig & Wolfe 1960, "Decomposition principle for linear programs"; finite convergence because there are finitely many extreme points and rays).
Bound hierarchy and the Lagrangian connection
For a MIP with integer blocks, convexification replaces each block's LP relaxation by $\operatorname{conv}(X_k)$, so
$$ z_{\mathrm{LP}} ;\le; z_{\mathrm{DW}} ;=; z_{\mathrm{LD}} ;\le; z_{\mathrm{IP}}, $$
where $z_{\mathrm{LD}}$ is the Lagrangian dual obtained by dualizing exactly the linking constraints. The first inequality is strict precisely when some $\operatorname{conv}(X_k)$ is strictly inside the block's LP relaxation — i.e., when the subproblems lack the integrality property (Geoffrion 1974, "Lagrangean relaxation for integer programming"). If every block prices integrally already (network flow, assignment, transportation blocks), then $z_{\mathrm{DW}} = z_{\mathrm{LP}}$ and DW buys structure and parallelism but no bound.
At any iteration with duals $\pi$, weak duality gives the Lagrangian bound
$$ L(\pi) ;=; \pi^{\top} b + \sum_{k} \min_{x \in X_k} (c_k - A_k^{\top}\pi)^{\top} x ;=; z_{\mathrm{RMP}} + \sum_{k} \bar{c}k ;\le; z{\mathrm{DW}}, $$
using LP duality on the RMP ($z_{\mathrm{RMP}} = \pi^{\top}b + \sum_k \sigma_k$). Track $\max_\pi L(\pi)$ across iterations: it sandwiches $z_{\mathrm{DW}}$ from below and supports early stopping long before pricing fully converges. If any block pricing is unbounded, $L(\pi) = -\infty$ for that $\pi$ and the bound update is skipped.
When to decompose — decision guidance
| Situation | Recommendation |
|---|---|
| Many blocks, few linking rows, integer blocks without integrality property | DW: bound improves and pricing decomposes |
| Blocks are pure network/assignment structures | DW bound = LP bound; decompose only for size/parallelism |
| Identical blocks (rolls, vehicles, machines) | DW with aggregated master; removes block symmetry |
| Linking rows dominate the matrix | Do not decompose; consider Benders if variables (not rows) split instead |
| Pricing has no exploitable structure and is itself a hard MIP | Expect pricing to dominate runtime; try heuristic pricing or skip DW |
| Compact LP already solves in seconds with a small gap | Skip DW; the overhead will not pay |
Detecting Block-Angular Structure
Before reformulating, confirm the structure actually decomposes. Remove the candidate linking rows and compute connected components of the row-column bipartite graph: each component is one block. The greedy density heuristic below is cheap and works when coupling rows are dense; production-grade detection uses hypergraph partitioning (Bergner et al. 2015, "Automatic Dantzig-Wolfe reformulation of mixed integer programs", the basis of GCG's detector).
"""Block-angular structure detection on the constraint matrix."""
import numpy as np
import scipy.sparse as sp
from scipy.sparse.csgraph import connected_components
def find_blocks(A: sp.csr_matrix,
linking_rows: np.ndarray) -> list[tuple[np.ndarray, np.ndarray]]:
"""Split A into independent blocks after removing the linking rows.
Returns one (row_indices, col_indices) pair per block, with row indices
referring to the original matrix. Columns appearing only in linking rows
form their own zero-row blocks and are reported too.
"""
m, n = A.shape
keep = np.setdiff1d(np.arange(m), np.asarray(linking_rows, dtype=int))
pattern = (A[keep, :] != 0).astype(np.int8)
# Bipartite graph: kept rows are nodes 0..len(keep)-1, columns follow.
graph = sp.bmat([[None, pattern], [pattern.T, None]], format="csr")
_, labels = connected_components(graph, directed=False)
blocks = []
for comp in np.unique(labels[len(keep):]): # components holding a column
rows = keep[labels[: len(keep)] == comp]
cols = np.where(labels[len(keep):] == comp)[0]
blocks.append((rows, cols))
return blocks
def greedy_linking_rows(A: sp.csr_matrix, max_linking: int) -> np.ndarray:
"""Cheap detection heuristic: peel off the densest rows one at a time.
Dense rows are the most likely coupling rows. Keep the smallest removal
set that maximizes the block count within the budget. Serious detection
uses hypergraph partitioning (Bergner et al. 2015).
"""
density = np.diff(A.tocsr().indptr) # nonzeros per row
order = np.argsort(-density)
best = np.array([], dtype=int)
best_count = len(find_blocks(A, best))
for take in range(1, max_linking + 1):
cand = order[:take]
count = len(find_blocks(A, cand))
if count > best_count:
best, best_count = cand, count
return np.sort(best)
if __name__ == "__main__":
# 3 blocks of 2 vars each, 2 block rows per block, 2 dense linking rows on top.
block = np.array([[3.0, 5.0], [4.0, 2.0]])
A = sp.csr_matrix(np.vstack([np.ones((2, 6)),
sp.block_diag([block] * 3).toarray()]))
linking = greedy_linking_rows(A, max_linking=3)
print("linking rows:", linking)
for rows, cols in find_blocks(A, linking):
print("block rows", rows, "cols", cols)
# Expected: linking rows [0 1]; three blocks with rows [2 3], [4 5], [6 7]
# and columns [0 1], [2 3], [4 5] respectively.
If the heuristic finds no decomposition, the model may still decompose after assigning a small number of "ambiguous" rows to blocks by hand, or after a variable permutation suggested by domain knowledge (one block per machine, period, or scenario). When the matrix is born from a model you wrote, prefer annotating blocks at modeling time over rediscovering them numerically.
Generic Dantzig-Wolfe Implementation
The driver below is problem-independent. Each block supplies its cost vector, its linking-row coefficients, and a pricing oracle; the driver owns the RMP, the artificial columns, the dual extraction, the ray handling, and the Lagrangian bound.
DANTZIG-WOLFE COLUMN GENERATION (minimization)
Input: blocks k = 1..K with costs c_k, linking coefficients A_k, sets X_k;
linking RHS b with senses; tolerance tol.
1 build RMP: linking rows, one convexity row per block,
big-M artificial columns so the RMP starts feasible
2 repeat
3 solve RMP (LP) -> z_RMP, duals pi (linking), sigma_k (convexity)
4 for each block k:
5 solve pricing v_k = min { (c_k - A_k' pi)' x : x in X_k }
6 if unbounded with extreme ray r: rc = (c_k - A_k' pi)' r
7 else: rc = v_k - sigma_k
8 if rc < -tol: add the point/ray column to the RMP
9 if no block returned a ray:
10 best_LB = max(best_LB, z_RMP + sum_k rc_k) # Lagrangian bound
11 until no column was added or z_RMP - best_LB <= tol
12 z_DW = z_RMP; recover x_k = sum_p lambda_kp x_k^p + sum_q mu_kq r_k^q
13 fail loudly if any artificial is positive (original problem infeasible)
"""Generic Dantzig-Wolfe driver. Save as dw_framework.py for the worked examples."""
import math
from dataclasses import dataclass
from typing import Callable
import gurobipy as gp
import numpy as np
from gurobipy import GRB
PricingResult = tuple[float, np.ndarray, bool] # (value, point or ray, is_ray)
@dataclass
class Block:
"""Subproblem k: cost c_k, linking coefficients A_k, and a pricing oracle.
pricing(price) minimizes price @ x over X_k. It returns (value, x, False)
for an optimal extreme point, or (price @ r, r, True) for an extreme ray r
that proves the pricing problem is unbounded.
"""
c: np.ndarray # shape (n_k,): original block costs
A: np.ndarray # shape (m0, n_k): linking-row coefficients
pricing: Callable[[np.ndarray], PricingResult]
def make_block_pricing(D: np.ndarray, d: np.ndarray,
vtype: str = GRB.CONTINUOUS) -> Callable[[np.ndarray], PricingResult]:
"""Build a pricing oracle for X_k = {x >= 0 : D x <= d} (LP or integer)."""
model = gp.Model("pricing")
model.Params.OutputFlag = 0
model.Params.InfUnbdInfo = 1 # make extreme rays available
model.Params.DualReductions = 0 # distinguish INFEASIBLE from UNBOUNDED
x = model.addMVar(D.shape[1], lb=0.0, vtype=vtype, name="x")
model.addConstr(D @ x <= d, name="block")
def pricing(price: np.ndarray) -> PricingResult:
"""Minimize price @ x over the block polyhedron."""
model.setObjective(price @ x, GRB.MINIMIZE)
model.optimize()
if model.Status == GRB.OPTIMAL:
return model.ObjVal, x.X.copy(), False
if model.Status == GRB.UNBOUNDED and vtype == GRB.CONTINUOUS:
ray = np.array(model.getAttr("UnbdRay", model.getVars()))
ray /= np.abs(ray).max() # scale for numerical sanity
return float(price @ ray), ray, True
raise RuntimeError(f"pricing status {model.Status}: integer blocks must be bounded")
return pricing
class DantzigWolfe:
"""Dantzig-Wolfe master with one convexity constraint per block."""
def __init__(self, blocks: list[Block], b: np.ndarray, senses: str,
big_m: float = 1e6) -> None:
"""senses is one character per linking row, each of '<', '=', '>'."""
self.blocks = blocks
self.b = np.asarray(b, dtype=float)
self.master = gp.Model("rmp")
self.master.Params.OutputFlag = 0
smap = {"<": GRB.LESS_EQUAL, "=": GRB.EQUAL, ">": GRB.GREATER_EQUAL}
self.link = [self.master.addLConstr(gp.LinExpr(), smap[s], rhs, name=f"link[{i}]")
for i, (s, rhs) in enumerate(zip(senses, self.b))]
self.conv = [self.master.addLConstr(gp.LinExpr(), GRB.EQUAL, 1.0, name=f"conv[{k}]")
for k in range(len(blocks))]
# Artificial variables keep the RMP feasible before real columns exist.
self.artificials: list[gp.Var] = []
for i, s in enumerate(senses):
if s in ("=", ">"):
self.artificials.append(self.master.addVar(
obj=big_m, name=f"art_up[{i}]", column=gp.Column([1.0], [self.link[i]])))
if s in ("=", "<"):
self.artificials.append(self.master.addVar(
obj=big_m, name=f"art_dn[{i}]", column=gp.Column([-1.0], [self.link[i]])))
for k in range(len(blocks)):
self.artificials.append(self.master.addVar(
obj=big_m, name=f"art_conv[{k}]", column=gp.Column([1.0], [self.conv[k]])))
self.columns: list[tuple[int, np.ndarray, bool, gp.Var]] = []
self.history: list[dict] = []
def _add_column(self, k: int, x: np.ndarray, is_ray: bool) -> None:
"""Append one point or ray column for block k to the RMP."""
blk = self.blocks[k]
coeffs = [float(v) for v in blk.A @ x]
constrs: list[gp.Constr] = list(self.link)
if not is_ray:
coeffs.append(1.0) # convexity coefficient (points only)
constrs.append(self.conv[k])
var = self.master.addVar(obj=float(blk.c @ x), lb=0.0,
name=f"col[{k},{len(self.columns)}]",
column=gp.Column(coeffs, constrs))
self.columns.append((k, x.copy(), is_ray, var))
def solve(self, tol: float = 1e-6, max_iters: int = 200) -> dict:
"""Run column generation; return bounds, recovered x_k, and the iteration log."""
best_lb = -math.inf
z_rmp = math.inf
for it in range(1, max_iters + 1):
self.master.optimize()
if self.master.Status != GRB.OPTIMAL:
raise RuntimeError(f"RMP status {self.master.Status}")
z_rmp = self.master.ObjVal
pi = np.array([c.Pi for c in self.link])
sigma = np.array([c.Pi for c in self.conv])
new_cols, rc_sum, ray_seen = 0, 0.0, False
for k, blk in enumerate(self.blocks):
value, x, is_ray = blk.pricing(blk.c - blk.A.T @ pi)
rc = value if is_ray else value - sigma[k]
ray_seen = ray_seen or is_ray
rc_sum += 0.0 if is_ray else rc
if rc < -tol:
self._add_column(k, x, is_ray)
new_cols += 1
if not ray_seen: # Lagrangian bound L(pi), see anatomy
best_lb = max(best_lb, z_rmp + rc_sum)
self.history.append({"iter": it, "z_rmp": z_rmp,
"lower_bound": best_lb, "new_cols": new_cols})
if new_cols == 0 or z_rmp - best_lb <= tol * max(1.0, abs(z_rmp)):
break
if any(a.X > 1e-6 for a in self.artificials):
raise RuntimeError("artificials positive at the end: instance infeasible "
"or big_m too small")
x_blocks = [np.zeros(len(blk.c)) for blk in self.blocks]
for k, x, _, var in self.columns:
x_blocks[k] += var.X * x
return {"z_dw": z_rmp, "lower_bound": best_lb,
"iterations": len(self.history), "x_blocks": x_blocks}
Design notes. The pricing oracle is a closure over a persistent Gurobi model, so each call only swaps the objective — no model rebuild. Reduced-cost sums feed the Lagrangian bound $L(\pi) = z_{\mathrm{RMP}} + \sum_k \bar{c}_k$ every iteration, which is what allows the early-exit gap test on line z_rmp - best_lb <= tol * .... Rays carry no convexity coefficient, exactly mirroring the reformulation. The recovered x_blocks solve the original constraints but are generally fractional in the original variables even when every generated column is integral — that is the convexification at work, and it is why integer answers need branch-and-price.
Worked Example 1: Multi-Plant Production Planning
Three plants produce two products in integer batches. Each plant k has its own two resource constraints (the block $X_k = {x \in \mathbb{Z}_+^2 : R_k x \le \mathrm{cap}_k}$), and the only coupling is demand coverage:
$$ \min \sum_{k=1}^{3}\sum_{j=1}^{2} c_{kj} x_{kj} \quad \text{s.t.} \quad \sum_{k} x_{kj} \ge \mathrm{dem}_j ;; (j=1,2), \qquad x_k \in X_k ;; (k=1,2,3). $$
Each plant's linking coefficients are $A_k = I_2$: plant k contributes its own production to each demand row. The blocks are bounded integer programs without the integrality property, so $z_{\mathrm{DW}}$ should strictly beat $z_{\mathrm{LP}}$.
"""Multi-plant production planning: compact MIP and its LP relaxation.
Save as multiplant.py; the DW run in the next block imports make_instance.
"""
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def make_instance() -> dict:
"""Three plants, two products, two resources per plant. Integer batches."""
return {
"cost": np.array([[4.0, 7.0], [5.0, 6.0], [6.0, 5.0]]), # (K, J)
"use": [np.array([[3.0, 5.0], [4.0, 2.0]]), # R_k: (2, J)
np.array([[2.0, 4.0], [5.0, 3.0]]),
np.array([[4.0, 3.0], [3.0, 4.0]])],
"cap": [np.array([25.0, 22.0]),
np.array([20.0, 24.0]),
np.array([26.0, 23.0])],
"demand": np.array([10.0, 8.0]),
}
def solve_compact(data: dict, relax: bool) -> float:
"""Solve the compact model; relax=True drops integrality on the batches."""
K, J = data["cost"].shape
model = gp.Model("multiplant")
model.Params.OutputFlag = 0
vtype = GRB.CONTINUOUS if relax else GRB.INTEGER
x = model.addMVar((K, J), lb=0.0, vtype=vtype, name="x")
for k in range(K):
model.addConstr(data["use"][k] @ x[k, :] <= data["cap"][k], name=f"capacity[{k}]")
model.addConstr(x.sum(axis=0) >= data["demand"], name="demand")
model.setObjective((data["cost"] * x).sum(), GRB.MINIMIZE)
model.optimize()
assert model.Status == GRB.OPTIMAL
return model.ObjVal
if __name__ == "__main__":
data = make_instance()
print(f"compact LP relaxation: {solve_compact(data, relax=True):.4f}")
print(f"compact IP optimum: {solve_compact(data, relax=False):.4f}")
# Expected: compact LP relaxation 90.0909, compact IP optimum 92.0000
Now decompose: one Block per plant with integer pricing over its own capacity polytope, identity linking coefficients, and > senses on the two demand rows.
"""Dantzig-Wolfe on the multi-plant instance: convexified integer blocks."""
import numpy as np
from gurobipy import GRB
from dw_framework import Block, DantzigWolfe, make_block_pricing # framework above
from multiplant import make_instance, solve_compact # compact model above
def run() -> None:
"""Compare z_LP <= z_DW <= z_IP on the three-plant instance."""
data = make_instance()
K, J = data["cost"].shape
blocks = [Block(c=data["cost"][k],
A=np.eye(J), # plant k adds x_kj to demand row j
pricing=make_block_pricing(data["use"][k], data["cap"][k],
vtype=GRB.INTEGER))
for k in range(K)]
dw = DantzigWolfe(blocks, b=data["demand"], senses=">" * J)
result = dw.solve()
z_lp = solve_compact(data, relax=True)
z_ip = solve_compact(data, relax=False)
print(f"z_LP = {z_lp:.4f} z_DW = {result['z_dw']:.4f} z_IP = {z_ip:.4f}")
print(f"iterations = {result['iterations']}, columns = {len(dw.columns)}")
for k, xk in enumerate(result["x_blocks"]):
print(f"plant {k}: x = {np.round(xk, 4)}")
for row in dw.history:
print(row)
if __name__ == "__main__":
run()
# Expected: z_LP = 90.0909 < z_DW = 91.4000 < z_IP = 92.0000
# (6 iterations, 12 columns; DW closes about 69% of the LP-IP gap).
# Recovered plan: plant 0 -> [5, 1], plant 1 -> [3.2, 2.6], plant 2 -> [1.8, 4.4];
# plants 1 and 2 are fractional convex combinations of integer pricing points.
Reading the run: the first two RMP values are dominated by big-M artificials; once real columns cover the demand rows, z_RMP drops to the 90s while the Lagrangian bound climbs (73.2 → 86.0 → 91.0 → 91.4) until the two meet at the DW optimum. Every column generated is an integer production plan for one plant, yet the master mixes them fractionally — the bound improved from 90.09 to 91.4 exactly because conv(X_k) cuts off fractional plans like x = (8.33, 0) that the compact LP relaxation allows.
Worked Example 2: Cutting Stock — Compact vs. DW Reformulation
The one-dimensional cutting-stock problem is the canonical demonstration that the same problem under two formulations gives very different bounds. The compact (Kantorovich) model has one identical block per stock roll $k$: $y_k \in {0,1}$ (roll used) and $x_{ik} \in \mathbb{Z}+$ (copies of item i cut from roll k), with block constraint $\sum_i l_i x{ik} \le W y_k$ and linking constraints $\sum_k x_{ik} \ge d_i$. Its LP relaxation is famously weak — exactly the material bound $\sum_i l_i d_i / W$ (de Carvalho 2002, "LP models for bin packing and cutting stock problems").
Applying DW with one convexity row per roll and then aggregating the identical blocks (sum the K convexity constraints into $\sum_p \lambda_p \le K$, which is slack at the optimum and can be dropped because the objective already counts used rolls) yields the Gilmore-Gomory pattern formulation (Gilmore & Gomory 1961, "A linear programming approach to the cutting-stock problem"):
$$ z_{\mathrm{DW}} = \min \sum_{p} \lambda_p \quad \text{s.t.} \quad \sum_{p} a_{ip}, \lambda_p \ge d_i ;;[\pi_i], \qquad \lambda \ge 0, $$
with one column per feasible cutting pattern $a \in \mathbb{Z}+^{n}$, $\sum_i l_i a_i \le W$. Pricing is an unbounded knapsack: a pattern improves the master iff $\sum_i \pi_i a_i > 1$. The DW bound is so strong that the modified integer round-up property $z{\mathrm{IP}} \le \lceil z_{\mathrm{DW}} \rceil + 1$ is conjectured to always hold (Scheithauer & Terno 1995).
"""Cutting stock, compact (Kantorovich) model. Save as kantorovich.py."""
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def make_instance() -> dict:
"""Tiny cutting-stock instance: roll width 10, three item lengths."""
return {"W": 10.0,
"lengths": np.array([6.0, 5.0, 3.0]),
"demand": np.array([2, 2, 4])}
def solve_kantorovich(data: dict, relax: bool) -> float:
"""One block per roll: y_k roll used, x_ik copies of item i cut from roll k."""
lengths, demand, W = data["lengths"], data["demand"], data["W"]
n_items, n_rolls = len(lengths), int(demand.sum()) # demand.sum() rolls suffice
model = gp.Model("kantorovich")
model.Params.OutputFlag = 0
if relax:
y = model.addMVar(n_rolls, lb=0.0, ub=1.0, name="y")
x = model.addMVar((n_items, n_rolls), lb=0.0, name="x")
else:
y = model.addMVar(n_rolls, vtype=GRB.BINARY, name="y")
x = model.addMVar((n_items, n_rolls), lb=0.0, vtype=GRB.INTEGER, name="x")
for k in range(n_rolls):
model.addConstr(lengths @ x[:, k] <= W * y[k], name=f"width[{k}]") # block k
model.addConstr(x.sum(axis=1) >= demand, name="demand") # linking
model.setObjective(y.sum(), GRB.MINIMIZE)
model.optimize()
assert model.Status == GRB.OPTIMAL
return model.ObjVal
if __name__ == "__main__":
data = make_instance()
print(f"Kantorovich LP bound: {solve_kantorovich(data, relax=True):.4f}")
print(f"material bound: {data['lengths'] @ data['demand'] / data['W']:.4f}")
print(f"IP optimum: {solve_kantorovich(data, relax=False):.4f}")
# Expected: Kantorovich LP bound 3.4000 (= material bound 3.4000), IP optimum 4.0000
The DW reformulation, with dynamic-programming knapsack pricing instead of a MIP oracle — this is the structure-exploitation payoff the Initial Assessment asks about:
"""Cutting stock, Dantzig-Wolfe (Gilmore-Gomory) reformulation with aggregation."""
import gurobipy as gp
import numpy as np
from gurobipy import GRB
from kantorovich import make_instance, solve_kantorovich # compact model above
def knapsack_pricing(lengths: np.ndarray, W: float, pi: np.ndarray) -> tuple[float, np.ndarray]:
"""Unbounded knapsack DP: max pi @ a s.t. lengths @ a <= W, a integer >= 0."""
cap = int(W)
sizes = lengths.astype(int)
value = np.zeros(cap + 1)
choice = np.full(cap + 1, -1, dtype=int)
for w in range(1, cap + 1):
value[w] = value[w - 1] # leave one unit of width unused
for i, s in enumerate(sizes):
if s <= w and value[w - s] + pi[i] > value[w]:
value[w] = value[w - s] + pi[i]
choice[w] = i
pattern = np.zeros(len(lengths))
w = cap
while w > 0:
if choice[w] >= 0:
pattern[choice[w]] += 1
w -= sizes[choice[w]]
else:
w -= 1 # unused width unit
return float(value[cap]), pattern
def solve_gilmore_gomory(data: dict, tol: float = 1e-9) -> dict:
"""Aggregated DW master: min sum(lambda), pattern coverage >= demand."""
lengths, demand, W = data["lengths"], data["demand"], data["W"]
n_items = len(lengths)
master = gp.Model("gg_master")
master.Params.OutputFlag = 0
cover = [master.addLConstr(gp.LinExpr(), GRB.GREATER_EQUAL, float(demand[i]),
name=f"cover[{i}]") for i in range(n_items)]
patterns: list[np.ndarray] = []
lam: list[gp.Var] = []
def add_pattern(a: np.ndarray) -> None:
"""Register pattern a as a master column with unit cost."""
col = gp.Column([float(a[i]) for i in range(n_items)], cover)
lam.append(master.addVar(obj=1.0, lb=0.0, name=f"lam[{len(patterns)}]", column=col))
patterns.append(a.copy())
for i in range(n_items): # initial single-item patterns
a = np.zeros(n_items)
a[i] = np.floor(W / lengths[i])
add_pattern(a)
iterations = 0
while True:
iterations += 1
master.optimize()
assert master.Status == GRB.OPTIMAL
pi = np.array([c.Pi for c in cover])
best_value, pattern = knapsack_pricing(lengths, W, pi)
if best_value <= 1.0 + tol: # reduced cost 1 - best >= 0
break
add_pattern(pattern)
z_dw = master.ObjVal
for v in lam: # integer solve on final patterns
v.VType = GRB.INTEGER
master.optimize()
z_int = master.ObjVal
used = [(patterns[j], round(lam[j].X)) for j in range(len(lam)) if lam[j].X > 0.5]
return {"z_dw": z_dw, "z_int_over_patterns": z_int,
"iterations": iterations, "n_patterns": len(patterns), "used": used}
if __name__ == "__main__":
data = make_instance()
res = solve_gilmore_gomory(data)
print(f"Kantorovich LP bound: {solve_kantorovich(data, relax=True):.4f}")
print(f"DW (pattern) LP bound: {res['z_dw']:.4f}")
print(f"integer over final patterns: {res['z_int_over_patterns']:.0f}")
print(f"Kantorovich IP optimum: {solve_kantorovich(data, relax=False):.0f}")
print(f"CG iterations: {res['iterations']}, patterns: {res['n_patterns']}")
for a, count in res["used"]:
print(f" cut {count} roll(s) with pattern {a.astype(int)}")
# Expected: Kantorovich LP 3.4000 < DW pattern LP 3.6667; integer solve over
# the final patterns gives 4 = IP optimum, in 2 CG iterations with 4 patterns.
# Used patterns: 1 x [0 2 0], 1 x [0 0 3], 2 x [1 0 1].
The reformulation comparison in one line: identical bounded integer blocks, dualized demand rows, knapsack pricing — and the bound jumps from 3.4 to 3.667 while the model shrinks from n_items × n_rolls integer variables to a handful of pattern columns. The integer solve over the generated columns is a heuristic (columns optimal for the LP need not contain an optimal integer solution); here it happens to hit the optimum, which the Kantorovich IP confirms. Note also the relation to the modified round-up property: $\lceil 3.667 \rceil = 4 = z_{\mathrm{IP}}$.
Advanced Techniques
Dual stabilization with a soft box
Early RMP duals are wild — extreme points of a dual polyhedron that barely constrains them — so pricing chases noise ("heading-in") and convergence crawls near the end ("tailing-off"). A soft box (du Merle et al. 1999, "Stabilized column generation"; the hard-box ancestor is Marsten, Hogan & Blankenship 1975, the Boxstep method) penalizes duals for leaving a trust region around a center, typically the best Lagrangian point found so far:
"""Soft dual box (du Merle et al. 1999) for stabilizing the DW master."""
import gurobipy as gp
import numpy as np
from gurobipy import GRB
def add_dual_box(master: gp.Model, link_constrs: list[gp.Constr],
center: np.ndarray, half_width: float,
eps: float) -> list[gp.Var]:
"""Pull each linking dual toward [center_i - w, center_i + w] via slack columns.
In the RMP dual, the +1 column with cost center_i + w enforces
pi_i <= center_i + w while the slack stays inside its bound; the -1 column
enforces pi_i >= center_i - w. The upper bound eps makes the box soft:
duals may leave it once a slack hits eps, so the master cannot become
unbounded when the box misses the true dual region. Shrink w and recenter
on the best Lagrangian point as CG converges; remove the columns and
re-optimize before declaring the master optimal.
"""
added: list[gp.Var] = []
for i, con in enumerate(link_constrs):
up = master.addVar(obj=float(center[i]) + half_width, lb=0.0, ub=eps,
name=f"stab_up[{i}]", column=gp.Column([1.0], [con]))
dn = master.addVar(obj=-(float(center[i]) - half_width), lb=0.0, ub=eps,
name=f"stab_dn[{i}]", column=gp.Column([-1.0], [con]))
added.extend([up, dn])
return added
def remove_stabilization(master: gp.Model, stab_vars: list[gp.Var]) -> None:
"""Drop the box columns; re-optimize before trusting duals or the bound."""
for v in stab_vars:
master.remove(v)
master.update()
if __name__ == "__main__":
m = gp.Model("toy_rmp")
m.Params.OutputFlag = 0
x = m.addVar(obj=5.0, name="x")
c = m.addLConstr(2.0 * x, GRB.GREATER_EQUAL, 4.0, name="link")
stab = add_dual_box(m, [c], center=np.array([1.0]), half_width=1.0, eps=10.0)
m.optimize()
print(f"dual with box: {c.Pi:.1f}")
remove_stabilization(m, stab)
m.optimize()
print(f"dual without box: {c.Pi:.1f}")
# Expected: dual with box 2.0 (capped at center + half_width), without box 2.5
A cheaper alternative that often suffices: solve the RMP with the barrier method and no crossover (Method=2, Crossover=0), which returns well-centered interior duals instead of vertex duals. Smoothing schemes (Wentges 1997; Pessoa et al. 2018 give the modern automatic variant) average the current duals with the best Lagrangian point and pair naturally with the bound tracking already in the framework. Full treatment belongs to the column-generation skill.
Convexification versus discretization
Everything above is convexification: $\lambda$ continuous, integrality conceptually retained on the original $x_k$. The alternative, discretization (Vanderbeck 2000; Vanderbeck & Savelsbergh 2006, "A generic view of Dantzig-Wolfe decomposition in mixed integer programming"), enumerates the integer points of $X_k$ (not just extreme points) and requires $\lambda$ integral in the master. Both give the same LP bound, but discretization matters when you must express integrality in master variables: with aggregated identical blocks, "$\lambda_p$ = number of blocks using pattern p" is meaningful and branchable — exactly the Gilmore-Gomory integer master — whereas convexified $\lambda$ has no integer meaning per block. Rule of practice: pure-integer bounded blocks plus aggregation calls for discretization semantics; mixed-integer or unbounded blocks stay with convexification.
Aggregating identical blocks
When all K blocks share the same data, per-block convexity constraints create a master in which permuting block indices yields equivalent solutions — heavy degeneracy and symmetry. Replace $\sum_p \lambda_{kp} = 1$ for each k by one aggregated constraint $\sum_p \lambda_p = K$ (or $\le K$, droppable when the objective counts active blocks, as in cutting stock). The aggregated master is smaller, prices once per iteration instead of K times, and removes the symmetry. To recover a per-block solution afterwards, greedily split the aggregated $\lambda$ into K unit-sum groups; any split is valid because the blocks are interchangeable.
Choosing DW-CG versus subgradient on the same dual
DW column generation and Lagrangian relaxation optimize the same dual function $L(\pi)$ — DW by an LP master holding all observed cuts of the dual function, subgradient methods by first-order steps. DW gives exact dual optima, a primal (convexified) solution for free, and clean stopping criteria, but each iteration costs an LP solve and master degeneracy can hurt. Subgradient iterations are nearly free and memory-light but need step-size tuning and never certify optimality tightly. A robust hybrid: run a few hundred subgradient steps to get a good dual point cheaply, then start DW-CG with the box of the section above centered there. Cross-check: at termination both must report the same bound up to tolerance — a mismatch means a bug in pricing or in the dualization.
Practical Challenges
The first RMP iterations return absurd objective values and duals. That is the big-M artificial phase: with no real columns, artificials carry the constraints and the duals equal ±big_m. It is harmless if big_m exceeds any realistic dual, but a too-small big_m silently produces an infeasible "optimum". Safer variants: a phase-1 master that minimizes artificial mass only, or seeding each block with one column priced at $\pi = 0$ (i.e., the block solved with original costs) so convexity rows are covered from the start.
Column generation stalls near the optimum (tailing-off). Dozens of iterations gain 0.001 each. Use the Lagrangian bound already computed by the framework: stop when $z_{\mathrm{RMP}} - \max L(\pi)$ falls below the accuracy you actually need — for a MIP bound, below 1 ulp of the rounding unit (e.g., gap < 1 when the objective is integral means the bound is proven). Stabilize the duals before blaming the method.
Your chosen "blocks" turn out not to be independent. A single forgotten coupling row makes pricing solutions invalid as columns. Always validate with find_blocks on the actual matrix: the union of returned blocks must cover every column, and every non-linking row must land in exactly one block. If a row straddles two would-be blocks, either promote it to linking (master grows) or merge the two blocks (pricing grows). Prefer merging when the merged pricing stays polynomial; prefer promoting when the row is one of few.
Integer pricing dominates the runtime. Each iteration solves K MIPs. Mitigations in order of payoff: replace the MIP oracle with a specialized algorithm (DP knapsack as in the cutting-stock example, labeling shortest path); accept any negative-reduced-cost heuristic column and call the exact oracle only when the heuristic fails (the final exact pass still certifies optimality); price blocks in parallel; keep a column pool and re-check old columns before re-pricing.
The DW bound equals the LP bound — all that reformulation for nothing. The blocks have the integrality property (Geoffrion 1974): e.g., assignment or pure network-flow blocks solve integrally as LPs. Re-draw the decomposition line so the hard integer structure sits inside the blocks: dualize a different constraint family, or enlarge blocks until they lose integrality. If no such split exists, DW still helps only via size/parallelism — say so explicitly rather than promising a better bound.
The master solution is fractional but the user needs an integer plan. Re-solving the final RMP with integer $\lambda$ (as in the cutting-stock example) is a quick heuristic but carries no guarantee — it can even be infeasible for set-partitioning masters. The exact route is branch-and-price with branching rules compatible with pricing; hand off to the column-generation skill rather than improvising branching on $\lambda$ variables, which destroys the pricing structure.
Identical blocks create a degenerate, symmetric master. K identical convexity rows admit factorially many equivalent bases; the simplex cycles through them and duals jump around. Aggregate the convexity constraints
…(truncated)