Extra material: Cutting stock#
The cutting stock problem is familiar to anyone who has cut parts out of stock materials. In the one-dimensional case, the stock materials are available in predetermined lengths and prices. The task is to cut a specific list of parts from the stock materials. The problem is to determine which parts to cut from each piece of stock material to minimize cost. This problem applies broadly to commercial applications, including the allocation of non-physical resources like capital budgeting or resource allocation.
This notebook presents several models and solution algorithms for the cutting stock problem.
Preamble: Install Pyomo and solvers#
We use the open-source solvers HiGHS for the linear models and SCIP for the bilinear mixed-integer models. On Colab, the following cell installs Pyomo, HiGHS and PySCIPOpt, which provides SCIP through its Python interface. For a local installation, run %pip install "pyomo>=6.10.1" highspy pyscipopt in a notebook first. The scip_direct interface avoids a separate solver executable and an IDAES download.
import sys
if "google.colab" in sys.modules:
%pip install -q "pyomo>=6.10.1" highspy pyscipopt
import pyomo.environ as pyo
MILO_solver = "appsi_highs"
MINLO_solver = "scip_direct"
MILO_SOLVER = pyo.SolverFactory(MILO_solver)
MINLO_SOLVER = pyo.SolverFactory(MINLO_solver)
assert MILO_SOLVER.available(), f"Solver {MILO_solver} is not available."
assert MINLO_SOLVER.available(), f"Solver {MINLO_solver} is not available."
Problem formulation#
Consider a set \({S}\) of available stock materials that can be cut to size to produce a set of finished parts. Each stock \(s\in S\) is characterized by a length \(l^S_s\), a cost \(c_s\) per piece, and is available in unlimited quantity. A customer order is received to product a set of finished products \(F\). Each finished product \(f\in F\) is specified by a required number \(d_f\) and length \(l^F_f\).
The cutting stock problem is to find a minimum cost solution to fulfill the customer order from the stock materials. The problem is illustrated is by data for an example given in the original paper by Gilmore and Gomory (1961).
Stocks
stocks |
length |
cost |
|---|---|---|
A |
5 |
6 |
B |
6 |
7 |
C |
9 |
10 |
Finished Parts
finished parts |
length |
demand |
|---|---|---|
S |
2 |
20 |
M |
3 |
10 |
L |
4 |
20 |
This information is represented in Python as nested dictionaries where the names for stocks and finished parts are used as indices.
stocks = {
"A": {"length": 5, "cost": 6},
"B": {"length": 6, "cost": 7},
"C": {"length": 9, "cost": 10},
}
finish = {
"S": {"length": 2, "demand": 20},
"M": {"length": 3, "demand": 10},
"L": {"length": 4, "demand": 20},
}
Patterns#
One approach to solving this problem is to create a list of all finished parts, a list of stocks for each length, and then use a set of binary decision variables to assign each finished product to a particular piece of stock. This approach will work well for a small problems, but the computational complexity scales much too rapidly with the size of the problem to be practical for business applications.
To address the issue of computational complexity, in 1961 Gilmore and Gomory introduced an additional data structure for the problem that is now referred to as “patterns”. A pattern is a list of finished parts that can be cut from a particular stock item.
A pattern \(p\) is specified by the stock \(s_p\) assigned to the pattern and integers \(a_{pf}\) that specify how many finished parts of type \(f\) are cut from stock \(s_p\). A pattern \(p\in P\) is feasible if
The function make_patterns defined below produces a partial list of feasible patterns for given sets of stocks and finished parts. Each pattern is represented as dictionary that specifies an associated stock item, and a dictionary of cuts that specify the finished parts cut from the stock. The algorithm is simple, it just considers every finished parts and stock items, then reports the number of parts \(f\) that can be cut from stock item \(s\).
def make_patterns(stocks, finish):
"""
Generates patterns of feasible cuts from stock lengths to meet specified finish lengths.
Parameters:
stocks (dict): A dictionary where keys are stock identifiers and values are dictionaries
with key 'length' representing the length of each stock.
finish (dict): A dictionary where keys are finish identifiers and values are dictionaries
with key 'length' representing the required finish lengths.
Returns:
patterns (list): A list of dictionaries, where each dictionary represents a pattern of cuts.
Each pattern dictionary contains 'stock' (the stock identifier) and 'cuts'
(a dictionary where keys are finish identifiers and the value is the number
of cuts from the stock for each finish).
"""
patterns = []
for f in finish:
feasible = False
for s in stocks:
# max number of f that fit on s
num_cuts = int(stocks[s]["length"] / finish[f]["length"])
# make pattern and add to list of patterns
if num_cuts > 0:
feasible = True
cuts_dict = {key: 0 for key in finish.keys()}
cuts_dict[f] = num_cuts
patterns.append({"stock": s, "cuts": cuts_dict})
if not feasible:
print(f"No feasible pattern was found for {f}")
return []
return patterns
patterns = make_patterns(stocks, finish)
display(patterns)
[{'stock': 'A', 'cuts': {'S': 2, 'M': 0, 'L': 0}},
{'stock': 'B', 'cuts': {'S': 3, 'M': 0, 'L': 0}},
{'stock': 'C', 'cuts': {'S': 4, 'M': 0, 'L': 0}},
{'stock': 'A', 'cuts': {'S': 0, 'M': 1, 'L': 0}},
{'stock': 'B', 'cuts': {'S': 0, 'M': 2, 'L': 0}},
{'stock': 'C', 'cuts': {'S': 0, 'M': 3, 'L': 0}},
{'stock': 'A', 'cuts': {'S': 0, 'M': 0, 'L': 1}},
{'stock': 'B', 'cuts': {'S': 0, 'M': 0, 'L': 1}},
{'stock': 'C', 'cuts': {'S': 0, 'M': 0, 'L': 2}}]
The function plot_patterns, defined below, displays a graphical depiction of the list of patterns.
import matplotlib.pyplot as plt
def plot_patterns(stocks, finish, patterns):
# set up figure parameters
lw = 0.6
cmap = plt.get_cmap("tab10")
colors = {f: cmap(k % 10) for k, f in enumerate(finish.keys())}
fig, ax = plt.subplots(1, 1, figsize=(8, 0.05 + 0.4 * len(patterns)))
for k, pattern in enumerate(patterns):
# get stock key/name
s = pattern["stock"]
# plot stock as a grey background
y_lo = (-k - lw / 2, -k - lw / 2)
y_hi = (-k + lw / 2, -k + lw / 2)
ax.fill_between((0, stocks[s]["length"]), y_lo, y_hi, color="k", alpha=0.1)
# overlay finished parts
xa = 0
for f, n in pattern["cuts"].items():
for j in range(n):
xb = xa + finish[f]["length"]
ax.fill_between((xa, xb), y_lo, y_hi, alpha=1.0, color=colors[f])
ax.plot((xb, xb), (y_lo[0], y_hi[0]), "w", lw=1, solid_capstyle="butt")
ax.text(
(xa + xb) / 2,
-k,
f,
ha="center",
va="center",
fontsize=6,
color="w",
weight="bold",
)
xa = xb
# clean up axes
ax.spines[["top", "right", "left", "bottom"]].set_visible(False)
ax.set_yticks(
range(0, -len(patterns), -1),
[pattern["stock"] for pattern in patterns],
fontsize=8,
)
plt.tight_layout()
return ax
ax = plot_patterns(stocks, finish, patterns)
Optimal cutting using known patterns#
Given a list of patterns, the optimization problem is to compute how many copies of each pattern should be cut to meet the demand for finished parts at minimum cost.
Let the index \(s_p\) denote the stock specified by pattern \(p\), and let \(x_p\) denote the number of copies of pattern \(p\) cut from stock \(s_p\). For a given list of patterns, the minimum cost optimization problem is a mixed integer linear optimization (MILO) subject to meeting demand constraints for each finished item.
The following cell is a Pyomo implementation of this optimization model.
# Given dictionaries of stocks and finished parts, and a list of patterns,
# find minimum choice of patterns to cut
def cut_patterns(stocks, finish, patterns):
m = pyo.ConcreteModel()
m.S = pyo.Set(initialize=list(stocks.keys()))
m.F = pyo.Set(initialize=list(finish.keys()))
m.P = pyo.RangeSet(0, len(patterns) - 1)
m.x = pyo.Var(m.P, domain=pyo.NonNegativeIntegers)
@m.Objective(sense=pyo.minimize)
def cost(m):
s = {p: patterns[p]["stock"] for p in m.P}
c = {p: stocks[s[p]]["cost"] for p in m.P}
return sum(c[p] * m.x[p] for p in m.P)
@m.Constraint(m.F)
def demand(m, f):
a = {p: patterns[p]["cuts"][f] for p in m.P}
return sum(a[p] * m.x[p] for p in m.P) >= finish[f]["demand"]
results = MILO_SOLVER.solve(m)
pyo.assert_optimal_termination(results)
return [round(m.x[p]()) for p in m.P], m.cost()
x, cost = cut_patterns(stocks, finish, patterns)
The following function plot_nonzero_patterns is wrapper for plot_patterns that removes unused patterns from graphic, shows the number of times each pattern is used, and adds cost to the title.
def plot_nonzero_patterns(stocks, finish, patterns, x, cost):
k = [j for j, _ in enumerate(x) if _ > 0.5]
ax = plot_patterns(stocks, finish, [patterns[j] for j in k])
ticks = [
f"{round(x[k])} x {pattern['stock']}"
for k, pattern in enumerate(patterns)
if x[k] > 0.5
]
ax.set_yticks(range(0, -len(k), -1), ticks, fontsize=8)
ax.set_title(f"Cost = {round(cost,2)}", fontsize=10)
return ax
ax = plot_nonzero_patterns(stocks, finish, patterns, x, cost)
Cutting Stock Problem: Bilinear reformulation#
The cut_patterns model requires a known list of cutting patterns. This works well if the patterns comprising an optimal solution to the problem are known. But since they are not initially known, an optimization model is needed that simultaneously solves for an optimal patterns and the cutting list.
Let binary variable \(b_{sp}\in\mathbb{Z}_2\) denote the assignment of stock \(s\) to pattern \(p\), and let \(P = 0, 1, \ldots, N_p-1\) index a list of patterns. For sufficiently large \(N_p\), an optimal solution to the stock cutting problem is given by the model
Since there is no ordering of the patterns, without loss of generality an additional constraint can be added to reduce the symmetries present in the problem.
This is a challenging optimization problem with a cost objective that is bilinear with respect to the the decision variables \(b_{sp}\) and \(x_p\), and a set of constraints for the demand of finished parts that are bilinear in the decision variables \(a_{fp}\) and \(x_p\). Because the constraints are a lower bound on a positive sum of bilinear terms, a simple substitution to create rotated quadratic cones fails to produce a convex program.
The following model is a direct translation of the bilinear optimization model into Pyomo. A solution is attempted using a mixed-integer nonlinear optimization (MINLO) solver.
Assume all lengths and costs are positive, demands are non-negative integers, and each requested part fits at least one stock. An optimal solution needs no more than \(U=\sum_f d_f\) stock pieces: any used piece that contributes nothing to meeting demand can be removed. We therefore bound each \(x_p\) and their sum by \(U\). The maximum individual demand is not such a bound. Choosing fewer than \(U\) pattern slots restricts the model; the two-slot example below is a small illustration, not a general choice of \(N_p\).
def bilinear_cut_stock(stocks, finish, Np=None):
if Np is None:
Np = len(finish)
m = pyo.ConcreteModel()
m.S = pyo.Set(initialize=list(stocks.keys()))
m.F = pyo.Set(initialize=list(finish.keys()))
m.P = pyo.RangeSet(0, Np - 1)
# sum of all finished parts
f_total_demand = sum(finish[f]["demand"] for f in m.F)
# define a with bounded values
def a_bounds(m, f, p):
return (0, max([int(stocks[s]["length"] / finish[f]["length"]) for s in m.S]))
m.a = pyo.Var(m.F, m.P, domain=pyo.NonNegativeIntegers, bounds=a_bounds)
m.b = pyo.Var(m.S, m.P, domain=pyo.Binary)
m.x = pyo.Var(m.P, domain=pyo.NonNegativeIntegers, bounds=(0, f_total_demand))
# minimize cost
@m.Objective(sense=pyo.minimize)
def cost(m):
c = {s: stocks[s]["cost"] for s in m.S}
return sum(c[s] * m.b[s, p] * m.x[p] for s in m.S for p in m.P)
# assign one stock to each pattern
@m.Constraint(m.P)
def assign_each_stock_to_pattern(m, p):
return sum(m.b[s, p] for s in m.S) == 1
# pattern feasibility
@m.Constraint(m.P)
def feasible_pattern(m, p):
return sum(m.a[f, p] * finish[f]["length"] for f in m.F) <= sum(
m.b[s, p] * stocks[s]["length"] for s in m.S
)
# demand constraints
@m.Constraint(m.F)
def demand(m, f):
return sum(m.a[f, p] * m.x[p] for p in m.P) >= finish[f]["demand"]
# order the patterns to reduce symmetries
@m.Constraint(m.P)
def order(m, p):
if p == 0:
return pyo.Constraint.Skip
return m.x[p] >= m.x[p - 1]
# upper bound on sum of all patterns
@m.Constraint()
def max_patterns(m):
return sum(m.x[p] for p in m.P) <= f_total_demand
results = MINLO_SOLVER.solve(m)
pyo.assert_optimal_termination(results)
cost = m.cost()
x = [round(m.x[p]()) for p in m.P]
patterns = []
for p in m.P:
a = {f: round(m.a[f, p]()) for f in m.F}
patterns.append({"stock": [s for s in m.S if m.b[s, p]() > 0.5][0], "cuts": a})
return patterns, x, cost
patterns, x, cost = bilinear_cut_stock(stocks, finish, 2)
plot_nonzero_patterns(stocks, finish, patterns, x, cost)
<Axes: title={'center': 'Cost = 170.0'}>
Pattern Generation: Bilinear Model#
From limited testing, the bilinear model for the cutting stock problem appears to work well for small data sets, but does not scale well for larger problem instances, at least for the solvers included in the testing. This shouldn’t be surprising given the non-convex nature of the problem, the exclusive use of integer and binary decision variables, and a high degree of symmetry in the model equations.
So rather than attempt to solve the full problem all at once, the following model assumes an initial list of patterns has been determined, perhaps using the make_patterns function defined above, then attempts to generate one more pattern that further reduces the objective function. The result remains a non-convex, bilinear optimization problem, but with fewer binary decision variables and at most one bilinear term in the objective and constraints.
The function generate_pattern_bilinear is a direct Pyomo implementation that uses a MINLO solver to create one additional feasible pattern that could be added to the list of known patterns.
Both the existing pattern counts \(x_p\) and the new count \(x'\) are integers bounded by \(U=\sum_f d_f\). Keeping these counts integer lets this step improve a cutting plan even after the linear relaxation has converged. With continuous counts, that improvement would already have been exhausted by complete linear-dual pricing. When an incumbent cost is supplied, the function returns None unless the additional pattern yields a strictly cheaper integer plan, within a numerical tolerance.
def generate_pattern_bilinear(stocks, finish, patterns, incumbent_cost=None, tol=1e-6):
m = pyo.ConcreteModel()
m.S = pyo.Set(initialize=list(stocks.keys()))
m.F = pyo.Set(initialize=list(finish.keys()))
m.P = pyo.RangeSet(0, len(patterns) - 1)
total_demand = sum(finish[f]["demand"] for f in m.F)
m.x = pyo.Var(m.P, domain=pyo.NonNegativeIntegers, bounds=(0, total_demand))
m.xp = pyo.Var(domain=pyo.NonNegativeIntegers, bounds=(0, total_demand))
def a_bounds(m, f):
return (0, max([int(stocks[s]["length"] / finish[f]["length"]) for s in m.S]))
m.ap = pyo.Var(m.F, domain=pyo.NonNegativeIntegers, bounds=a_bounds)
m.bp = pyo.Var(m.S, domain=pyo.Binary)
@m.Objective(sense=pyo.minimize)
def cost(m):
s = {p: patterns[p]["stock"] for p in m.P}
c = {p: stocks[s[p]]["cost"] for p in m.P}
return sum(c[p] * m.x[p] for p in m.P) + m.xp * sum(
m.bp[s] * stocks[s]["cost"] for s in m.S
)
@m.Constraint()
def sum_to_one(m):
return sum(m.bp[s] for s in m.S) == 1
@m.Constraint()
def add_pattern(m):
return sum(m.ap[f] * finish[f]["length"] for f in m.F) <= sum(
m.bp[s] * stocks[s]["length"] for s in m.S
)
@m.Constraint(m.F)
def demand(m, f):
return (
sum(patterns[p]["cuts"][f] * m.x[p] for p in m.P) + m.ap[f] * m.xp
>= finish[f]["demand"]
)
results = MINLO_SOLVER.solve(m)
pyo.assert_optimal_termination(results)
if m.xp() < 0.5 or (incumbent_cost is not None and m.cost() >= incumbent_cost - tol):
return None
new_pattern = {
"stock": [s for s in m.S if m.bp[s]() > 0.5][0],
"cuts": {f: round(m.ap[f]()) for f in m.F},
}
return new_pattern
stocks = {
"A": {"length": 5, "cost": 6},
"B": {"length": 6, "cost": 7},
"C": {"length": 9, "cost": 10},
}
finish = {
"S": {"length": 2, "demand": 20},
"M": {"length": 3, "demand": 10},
"L": {"length": 4, "demand": 20},
}
patterns = make_patterns(stocks, finish)
generate_pattern_bilinear(stocks, finish, patterns)
{'stock': 'B', 'cuts': {'S': 1, 'M': 0, 'L': 1}}
Pattern Generation: Linear Dual#
A common approach to pattern generation for stock cutting begins by relaxing the MILO optimization problem with known patterns. The integer variables \(x_p\) are relaxed to non-negative reals.
Let \(\pi_f \geq 0\) be the dual variables associated with the demand constraints. A large positive value \(\pi_f\) suggests a high value for including finished part \(f\) in a new pattern. This motivates a set of dual optimization problems where the objective is to construct a new patterns that maximizes the the marginal value of each stock \(s\in S\).
The pattern demonstrating the largest return \(V_s\) is returned as a candidate to add the list of patterns.
Pricing stops when no stock has \(V_s>0\), within tolerance. The function returns the improving pattern (or None) and a lower bound from the relaxation. To keep the bound valid when pricing stops at a small positive tolerance, the dual prices are scaled down if needed so that no feasible pattern is valued above its stock cost.
import pyomo.environ as pyo
def generate_pattern_dual(stocks, finish, patterns, tol=1e-6):
m = pyo.ConcreteModel()
m.S = pyo.Set(initialize=list(stocks.keys()))
m.F = pyo.Set(initialize=list(finish.keys()))
m.P = pyo.RangeSet(0, len(patterns) - 1)
m.x = pyo.Var(m.P, domain=pyo.NonNegativeReals)
@m.Objective(sense=pyo.minimize)
def cost(m):
s = {p: patterns[p]["stock"] for p in m.P}
c = {p: stocks[s[p]]["cost"] for p in m.P}
return sum(c[p] * m.x[p] for p in m.P)
@m.Constraint(m.F)
def demand(m, f):
a = {p: patterns[p]["cuts"][f] for p in m.P}
return sum(a[p] * m.x[p] for p in m.P) >= finish[f]["demand"]
# solve and get duals
m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
results = MILO_SOLVER.solve(m)
pyo.assert_optimal_termination(results)
# solve for reduced cost for each stock
@m.Block(m.S)
def subprob(b, s):
m = b.model()
b.a = pyo.Var(m.F, domain=pyo.NonNegativeIntegers)
@b.Objective(sense=pyo.maximize)
def marginal_value(b):
return sum(b.a[f] * m.dual[m.demand[f]] for f in m.F)
@b.Constraint()
def stock_length(b):
return sum(b.a[f] * finish[f]["length"] for f in m.F) <= stocks[s]["length"]
results = MILO_SOLVER.solve(b)
pyo.assert_optimal_termination(results)
return b
reduced_value = [m.subprob[s].marginal_value() - stocks[s]["cost"] for s in m.S]
s = list(m.S)[reduced_value.index(max(reduced_value))]
scale = max(1.0, max(m.subprob[s].marginal_value() / stocks[s]["cost"] for s in m.S))
lower_bound = sum(finish[f]["demand"] * m.dual[m.demand[f]] for f in m.F) / scale
if max(reduced_value) <= tol:
return None, lower_bound
return {"stock": s, "cuts": {f: round(m.subprob[s].a[f]()) for f in m.F}}, lower_bound
stocks = {
"A": {"length": 5, "cost": 6},
"B": {"length": 6, "cost": 7},
"C": {"length": 9, "cost": 10},
}
finish = {
"S": {"length": 2, "demand": 20},
"M": {"length": 3, "demand": 10},
"L": {"length": 4, "demand": 20},
}
patterns = make_patterns(stocks, finish)
generate_pattern_dual(stocks, finish, patterns)
({'stock': 'C', 'cuts': {'S': 1, 'M': 1, 'L': 1}}, 168.74999999999997)
The following cell times one call to each method on a larger data set. These timings illustrate the work involved on the current machine; they are not a solver benchmark, and the two methods solve different subproblems.
stocks = {
'log': {'length': 100, 'cost': 1},
}
finish = {
1: {'length': 75.0, 'demand': 38},
2: {'length': 75.0, 'demand': 44},
3: {'length': 75.0, 'demand': 30},
4: {'length': 75.0, 'demand': 41},
5: {'length': 75.0, 'demand': 36},
6: {'length': 53.8, 'demand': 33},
7: {'length': 53.0, 'demand': 36},
8: {'length': 51.0, 'demand': 41},
9: {'length': 50.2, 'demand': 35},
10: {'length': 32.2, 'demand': 37},
11: {'length': 30.8, 'demand': 44},
12: {'length': 29.8, 'demand': 49},
13: {'length': 20.1, 'demand': 37},
14: {'length': 16.2, 'demand': 36},
15: {'length': 14.5, 'demand': 42},
16: {'length': 11.0, 'demand': 33},
17: {'length': 8.6, 'demand': 47},
18: {'length': 8.2, 'demand': 35},
19: {'length': 6.6, 'demand': 49},
20: {'length': 5.1, 'demand': 42},
}
patterns = make_patterns(stocks, finish)
print("Testing generate_patterns_bilinear: ", end="")
%timeit -n 1 -r 1 generate_pattern_bilinear(stocks, finish, patterns)
print("Testing generate_patterns_dual: ", end="")
%timeit -n 1 -r 1 generate_pattern_dual(stocks, finish, patterns)
Testing generate_patterns_bilinear: 820 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
Testing generate_patterns_dual: 61.7 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
A hybrid solution algorithm using pattern generation#
We start with make_patterns and use linear-dual pricing until no improving column remains. Solving the integer model over these patterns gives a feasible cutting plan and an upper bound; the priced relaxation supplies a lower bound. Phase 2 searches for a single additional pattern that strictly reduces the integer cost, repeating while an improvement is found.
The final plan is optimal over the generated patterns. A gap between its cost and the lower bound can remain: failure to improve by adding one pattern does not rule out improvements requiring several new patterns together. The algorithm reports both bounds instead of claiming that every instance is solved globally.
def cut_stock(stocks, finish):
patterns = make_patterns(stocks, finish)
if not patterns:
raise ValueError("Every requested part must fit at least one stock.")
print("Phase 1 ", end=".")
while True:
new_pattern, lower_bound = generate_pattern_dual(stocks, finish, patterns)
if new_pattern is None:
break
if new_pattern in patterns:
raise RuntimeError("Pricing returned an existing improving pattern; check tolerances.")
patterns.append(new_pattern)
print(end=".")
x, cost = cut_patterns(stocks, finish, patterns)
print(f" Cost = {cost:.2f}; LP lower bound = {lower_bound:.6f}")
print("Phase 2 ", end=".")
while True:
new_pattern = generate_pattern_bilinear(stocks, finish, patterns, incumbent_cost=cost)
if new_pattern is None:
break
if new_pattern in patterns:
raise RuntimeError("An existing pattern cannot improve an optimal restricted model.")
patterns.append(new_pattern)
x, improved_cost = cut_patterns(stocks, finish, patterns)
assert improved_cost < cost - 1e-6
cost = improved_cost
print(end=".")
print(f" Cost = {cost:.2f}; bound gap = {max(0, cost - lower_bound):.6f}")
used = [i for i, count in enumerate(x) if count > 0]
return [patterns[i] for i in used], [x[i] for i in used], cost
stocks = {
"A": {"length": 5, "cost": 6},
"B": {"length": 6, "cost": 7},
"C": {"length": 9, "cost": 10},
}
finish = {
"S": {"length": 2, "demand": 20},
"M": {"length": 3, "demand": 10},
"L": {"length": 4, "demand": 20},
}
patterns, x, cost = cut_stock(stocks, finish)
plot_nonzero_patterns(stocks, finish, patterns, x, cost)
Phase 1 ... Cost = 170.00; LP lower bound = 170.000000
Phase 2 . Cost = 170.00; bound gap = 0.000000
<Axes: title={'center': 'Cost = 170.0'}>
Examples#
Example from JuMP documentation for column generation#
stocks = {
"log": {"length": 100, "cost": 1},
}
finish = {
1: {"length": 75.0, "demand": 38},
2: {"length": 75.0, "demand": 44},
3: {"length": 75.0, "demand": 30},
4: {"length": 75.0, "demand": 41},
5: {"length": 75.0, "demand": 36},
6: {"length": 53.8, "demand": 33},
7: {"length": 53.0, "demand": 36},
8: {"length": 51.0, "demand": 41},
9: {"length": 50.2, "demand": 35},
10: {"length": 32.2, "demand": 37},
11: {"length": 30.8, "demand": 44},
12: {"length": 29.8, "demand": 49},
13: {"length": 20.1, "demand": 37},
14: {"length": 16.2, "demand": 36},
15: {"length": 14.5, "demand": 42},
16: {"length": 11.0, "demand": 33},
17: {"length": 8.6, "demand": 47},
18: {"length": 8.2, "demand": 35},
19: {"length": 6.6, "demand": 49},
20: {"length": 5.1, "demand": 42},
}
patterns, x, cost = cut_stock(stocks, finish)
plot_nonzero_patterns(stocks, finish, patterns, x, cost)
Phase 1 ................ Cost = 334.00; LP lower bound = 334.000000
Phase 2 . Cost = 334.00; bound gap = 0.000000
<Axes: title={'center': 'Cost = 334.0'}>
Example from Wikipedia#
https://en.wikipedia.org/wiki/Cutting_stock_problem
The minimum number of rolls is 73.0.
stocks = {
"roll": {"length": 5600, "cost": 1},
}
finish = {
1380: {"length": 1380, "demand": 22},
1520: {"length": 1520, "demand": 25},
1560: {"length": 1560, "demand": 12},
1710: {"length": 1710, "demand": 14},
1820: {"length": 1820, "demand": 18},
1880: {"length": 1880, "demand": 18},
1930: {"length": 1930, "demand": 20},
2000: {"length": 2000, "demand": 10},
2050: {"length": 2050, "demand": 12},
2100: {"length": 2100, "demand": 14},
2140: {"length": 2140, "demand": 16},
2150: {"length": 2150, "demand": 18},
2200: {"length": 2200, "demand": 20},
}
patterns, x, cost = cut_stock(stocks, finish)
plot_nonzero_patterns(stocks, finish, patterns, x, cost)
Phase 1 ........................... Cost = 73.00; LP lower bound = 72.916667
Phase 2 . Cost = 73.00; bound gap = 0.083333
<Axes: title={'center': 'Cost = 73.0'}>
Woodworking: Problem data from Google sheets#
Find a minimum cost order of 2x4 lumber to build the “One Arm 2x4 Outdoor Sofa” described by Ana White.
Image source: www.ana-white.com
Data source: https://docs.google.com/spreadsheets/d/1ZX7KJ2kwTGgyqEv_a3LOG0nQSxsc38Ykk53A7vGWAFU/edit#gid=1104632299
This example reads a public sheet at run time; prices and orders may change, so a later run can differ from the saved results. The kerf setting is the material lost per cut. We charge one kerf per finished piece, including the last piece on a stock; this is a conservative convention that does not assume an offcut-free final cut.
import pandas as pd
def read_google_sheet(sheet_id, sheet_name):
"""
Reads a Google Sheet and returns a pandas DataFrame.
This function reads a Google Sheet with the specified sheet ID and sheet name,
and returns a pandas DataFrame with the data. The column names are converted to
lowercase.
Args:
sheet_id (str): The Google Sheet ID.
sheet_name (str): The name of the sheet to read.
Returns:
df (pd.DataFrame): A pandas DataFrame containing the data from the Google Sheet.
"""
url = f"https://docs.google.com/spreadsheets/d/{sheet_id}/gviz/tq?tqx=out:csv&sheet={sheet_name}"
df = pd.read_csv(url)
df.columns = map(str.lower, df.columns)
return df
# Google Sheet ID
sheet_id = "1ZX7KJ2kwTGgyqEv_a3LOG0nQSxsc38Ykk53A7vGWAFU"
# read settings
settings_df = read_google_sheet(sheet_id, "settings")
print("\nSettings")
display(settings_df)
# read parts
finish_df = read_google_sheet(sheet_id, "finish")
print("\nFinish")
display(finish_df)
# read and display stocks
stocks_df = read_google_sheet(sheet_id, "stocks")
if "price" not in stocks_df.columns:
raise ValueError("The stocks sheet must provide a price for each stock.")
print("\nStocks")
display(stocks_df)
Settings
| setting | value | |
|---|---|---|
| 0 | kerf | 0.125 |
Finish
| kind | length | quantity | label | |
|---|---|---|---|---|
| 0 | 2x4 | 70.50 | 3 | 70.50 |
| 1 | 2x4 | 25.50 | 10 | 25.50 |
| 2 | 2x4 | 12.50 | 1 | 12.50 |
| 3 | 2x4 | 72.00 | 6 | 72.00 |
| 4 | 2x4 | 70.75 | 1 | 70.75 |
| 5 | 2x4 | 28.50 | 1 | 28.50 |
Stocks
| kind | length | price | |
|---|---|---|---|
| 0 | 2x4 | 36 | 1.68 |
| 1 | 2x4 | 48 | 1.86 |
| 2 | 2x4 | 72 | 2.57 |
| 3 | 2x4 | 84 | 2.65 |
| 4 | 2x4 | 96 | 2.92 |
| 5 | 2x4 | 120 | 3.67 |
| 6 | 2x4 | 144 | 4.40 |
| 7 | 2x4 | 168 | 5.14 |
| 8 | 2x4 | 192 | 6.92 |
| 9 | 2x4 | 216 | 8.62 |
| 10 | 2x4 | 240 | 10.40 |
| 11 | 2x6 | 96 | 4.43 |
| 12 | 2x6 | 192 | 9.36 |
kinds = sorted(finish_df["kind"].unique())
kerf = float(settings_df.set_index("setting").loc["kerf", "value"])
purchase_rows = []
for kind in kinds:
print(f"Kind = {kind}")
finish = dict()
for i in finish_df.loc[finish_df["kind"] == kind].index:
finish[finish_df.loc[i, "label"]] = {
"length": finish_df.loc[i, "length"] + kerf,
"demand": finish_df.loc[i, "quantity"],
}
stocks = dict()
for i in stocks_df.loc[stocks_df["kind"] == kind].index:
stocks[stocks_df.loc[i, "length"]] = {
"length": stocks_df.loc[i, "length"],
"cost": stocks_df.loc[i, "price"],
}
patterns, x, cost = cut_stock(stocks, finish)
plot_nonzero_patterns(stocks, finish, patterns, x, cost)
for pattern, pieces in zip(patterns, x):
purchase_rows.append({"kind": kind, "stock length": pattern["stock"], "pieces": pieces})
Kind = 2x4
Phase 1 ...... Cost = 33.18; LP lower bound = 32.163333
Phase 2 .. Cost = 32.28; bound gap = 0.116667
Purchase list#
purchase_list = (
pd.DataFrame(purchase_rows)
.groupby(["kind", "stock length"], as_index=False)["pieces"].sum()
.sort_values(["kind", "stock length"])
)
purchase_list
| kind | stock length | pieces | |
|---|---|---|---|
| 0 | 2x4 | 144 | 5 |
| 1 | 2x4 | 168 | 2 |
References#
The one dimensional cutting stock problem addressed in this notebook is generally attributed to two classic papers by Gilmore and Gomory. This first paper considers the more general case of stocks available in multiple lengths, while the second paper specializes to the needs of a paper trimming operation.
Gilmore, P. C., & Gomory, R. E. (1961). A linear programming approach to the cutting-stock problem. Operations research, 9(6), 849-859. [jstor]
Gilmore, P. C., & Gomory, R. E. (1963). A linear programming approach to the cutting stock problem—Part II. Operations research, 11(6), 863-888. [jstor]
A useful survey of subsequent development of the cutting stock problem is given by:
Haessler, R. W., & Sweeney, P. E. (1991). Cutting stock problems and solution procedures. European Journal of Operational Research, 54(2), 141-150. [pdf]
Delorme, M., Iori, M., & Martello, S. (2016). Bin packing and cutting stock problems: Mathematical models and exact algorithms. European Journal of Operational Research, 255(1), 1-20. [sciencedirect]
The solution proposed by Gilmore and Gomory has been refined over time and now generally referred to as “column generation”. A number of tutorial implemenations are available, these are representative:
More recently, the essential bilinear structure of the problem has been noted, and various convex transformations of the problem have been studied:
Harjunkoski, I., Westerlund, T., Pörn, R., & Skrifvars, H. (1998). Different transformations for solving non-convex trim-loss problems by MINLP. European Journal of Operational Research, 105(3), 594-603. [abo.fi][sciencedirect]
Harjunkoski, I., Pörn, R., & Westerlund, T. (1999). Exploring the convex transformations for solving non-convex bilinear integer problems. Computers & Chemical Engineering, 23, S471-S474. [sciencedirect