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.

# install AMPL and solvers
%pip install -q amplpy

SOLVER_MILO = "highs"
SOLVER_MINLO = "ipopt"

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["coin", "highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 5.6/5.6 MB 17.3 MB/s eta 0:00:00
?25hUsing default Community Edition License for Colab. Get yours at: https://ampl.com/ce
Licensed to AMPL Community Edition License for the AMPL Model Colaboratory (https://colab.ampl.com).

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 Gamory (1961).

Stocks

stocks
\(s\)

length
\(l^S_s\)

cost
\(c_s\)

A

5

6

B

6

7

C

9

10

Finished Parts

finished parts
\(f\)

length
\(l^F_f\)

demand
\(d_f\)

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 Gamory 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

\[ \begin{align} \sum_{f\in F}a_{pf}l^F_f & \leq l^S_{s_p} \end{align} \]

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,
    )

    return ax


ax = plot_patterns(stocks, finish, patterns)
../../_images/a01b4cb77ecbe8898b5bdac444991c8abe20821028b88f01fce2fe28600cf27e.png

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_{s_p}\) denote the number pieces of stock \(s_p\) is used. 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.

\[\begin{split} \begin{align} \min\quad & \sum_{p\in P} c_{s_p} x_{s_p} \\ \text{s.t.}\quad & \sum_{p\in P}a_{pf} x_{s_p} \geq d_f && \forall f\in F\\ & x_{s_p} \in \mathbb{Z}_+ && \forall p\in P\\ \end{align} \end{split}\]

The following cell is an AMPL 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 = AMPL()

    m.eval(
        """
        set S;
        set F;
        set P;

        param c{P};
        param a{F, P};
        param demand_finish{F};

        var x{P} integer >= 0;

        minimize cost:
            sum{p in P} c[p]*x[p];

        subject to demand{f in F}:
            sum{p in P} a[f,p]*x[p] >= demand_finish[f];

    """
    )

    m.set["S"] = list(stocks.keys())
    m.set["F"] = list(finish.keys())
    m.set["P"] = list(range(len(patterns)))

    s = {p: patterns[p]["stock"] for p in range(len(patterns))}
    c = {p: stocks[s[p]]["cost"] for p in range(len(patterns))}
    m.param["c"] = c
    a = {
        (f, p): patterns[p]["cuts"][f]
        for p in range(len(patterns))
        for f in finish.keys()
    }
    m.param["a"] = a
    m.param["demand_finish"] = {
        f_part: finish[f_part]["demand"] for f_part in finish.keys()
    }

    m.option["solver"] = SOLVER_MILO
    m.get_output("solve;")

    return [m.var["x"][p].value() for p in range(len(patterns))], m.obj["cost"].value()


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]
    ax = plot_patterns(stocks, finish, [patterns[j] for j in k])
    ticks = [
        f"{x[k]} x {pattern['stock']}" for k, pattern in enumerate(patterns) if x[k] > 0
    ]
    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)
../../_images/ed5fcd8049eec227fb1e1add91eb20c9520f338de54d1f468fb99077c81097f4.png

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

\[\begin{split} \begin{align} \min\quad & \sum_{s\in S} \sum_{p\in P} c_{s} b_{sp} x_{p} \\ \text{s.t.}\quad & \sum_{s\in S}b_{sp} = 1 && \forall p\in P \\ & \sum_{f\in F}a_{fp}l^F_f \leq \sum_{s\in S} b_{sp} l^S_s && \forall p\in P \\ & \sum_{p\in P}a_{fp} x_{p} \geq d_f && \forall f\in F\\ & a_{fp}, x_p \in \mathbb{Z}_+ && \forall f\in F, \forall p\in P \\ & b_{sp} \in \mathbb{Z}_2 && \forall s\in S, \forall p\in P \\ \end{align} \end{split}\]

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.

\[ \begin{align} & x_{p-1} \geq x_{p} && \forall p\in P, p > 0 \end{align} \]

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 AMPL. A solution is attempted using a mixed-integer nonlinear optimization (MINLO) solver.

def bilinear_cut_stock(stocks, finish, Np=len(finish)):
    m = AMPL()

    m.eval(
        """
        set S;
        set F;
        set P;

        param c{S};

        # length stock
        param length_s{S};
        # length finished pieces
        param length_f{F};

        param demand_finish{F};
        param a_upper_bound{F};

        # sum of all finished parts
        param f_total_demand;

        var a{f in F, p in P} integer >= 0 <= a_upper_bound[f];
        var b{S, P} binary;
        var x{P} integer >= 0 <= f_total_demand;

        minimize cost:
            sum{p in P, s in S} c[s] * b[s,p] * x[p];

        subject to assign_each_stock_to_pattern{p in P}:
          sum{s in S} b[s,p] = 1;

        subject to feasible_pattern{p in P}:
          sum{f in F} a[f,p] * length_f[f] <= sum{s in S} b[s,p] * length_s[s];

        subject to demand{f in F}:
           sum{p in P} a[f,p]*x[p] >= demand_finish[f];

        # order the patterns to reduce symmetries
        subject to order{p in P : p >= 1}:
          x[p] <= x[p-1];

        subject to max_patterns:
          sum{p in P} x[p] <= f_total_demand;
    """
    )

    m.set["S"] = list(stocks.keys())
    m.set["F"] = list(finish.keys())
    m.set["P"] = list(range(Np))

    m.param["c"] = {s: stocks[s]["cost"] for s in list(stocks.keys())}
    m.param["length_s"] = {s: stocks[s]["length"] for s in stocks.keys()}
    m.param["length_f"] = {f: finish[f]["length"] for f in finish.keys()}
    m.param["demand_finish"] = {f: finish[f]["demand"] for f in finish.keys()}

    m.eval("let f_total_demand := max{f in F} demand_finish[f];")
    # or m.param['f_total_demand'] = max([finish[f]['demand'] for f in finish.keys()])
    a_upper_bound = {
        f: max([int(stocks[s]["length"] / finish[f]["length"]) for s in stocks.keys()])
        for f in finish.keys()
    }
    m.param["a_upper_bound"] = a_upper_bound

    m.option["solver"] = SOLVER_MINLO
    m.get_output("solve;")

    cost = m.obj["cost"].value()
    x = [round(m.var["x"][p].value()) for p in range(Np)]

    # Retrieve the patterns
    patterns = []
    for p in range(Np):
        a = {f: round(m.var["a"][f, p].value()) for f in finish.keys()}
        patterns.append(
            {
                "stock": [s for s in stocks.keys() if m.var["b"][s, p].value() > 0][0],
                "cuts": a,
            }
        )

    return patterns, x, cost


patterns, x, cost = bilinear_cut_stock(stocks, finish, 2)
plot_nonzero_patterns(stocks, finish, patterns, x, cost);
../../_images/40757b31f44cead31dd70e24ec184bb42e015740598fcd7e1e83cbbdab16b5a0.png

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.

\[\begin{split} \begin{align} \min\quad & \sum_{p\in P} c_{s_p} x_{s_p} + x' \sum_{s\in S} b_s c_s\\ \text{s.t.}\quad & \sum_{s\in S}b'_{s} = 1 \\ & \sum_{f\in F}a'_{f}l^F_f \leq \sum_{s\in S} b'_{s} l^S_s \\ & \sum_{p\in P}a_{fp} x_{s_p} + a'_f x'\geq d_f && \forall f\in F\\ & a'_{f}, x_p \in \mathbb{Z}_+ && \forall f\in F, \forall p\in P \\ & b'_{s} \in \mathbb{Z}_2 && \forall s\in S \\ \end{align} \end{split}\]

The function generate_pattern_bilinear is a direct AMPL implementation that uses a MINLO solver to create one additional feasible pattern that could be added to the list of known patterns.

pattern_bilinear_prob = AMPL()
pattern_bilinear_prob.option["solver"] = SOLVER_MINLO
pattern_bilinear_prob.eval(
    """
        set S;
        set F;
        set P;

        param c{P union S};

        # length stock
        param length_s{S};
        # length finished pieces
        param length_f{F};

        param demand_finish{F};
        param a{F, P};
        param ap_upper_bound{F};

        var ap{f in F} integer >= 0 <= ap_upper_bound[f];
        var bp{S} binary;
        var x{P} integer >= 0;
        var xp integer >= 0;

        minimize cost:
            sum{p in P} c[p] * x[p] + xp * sum{s in S} bp[s]*c[s];

        subject to assign_each_stock_to_pattern:
          sum{s in S} bp[s] = 1;

        subject to add_pattern:
          sum{f in F} ap[f] * length_f[f] <= sum{s in S} bp[s] * length_s[s];

        subject to demand{f in F}:
           sum{p in P} a[f,p]*x[p] + ap[f]*xp >= demand_finish[f];
"""
)


def generate_pattern_bilinear(stocks, finish, patterns):
    m = pattern_bilinear_prob
    m.eval("reset data;")

    m.set["S"] = list(stocks.keys())
    m.set["F"] = list(finish.keys())
    m.set["P"] = list(range(len(patterns)))

    s = {p: patterns[p]["stock"] for p in range(len(patterns))}
    c = {p: stocks[s[p]]["cost"] for p in range(len(patterns))}
    cstocks = {s: stocks[s]["cost"] for s in list(stocks.keys())}
    c.update(cstocks)
    m.param["c"] = c

    a = {
        (f, p): patterns[p]["cuts"][f]
        for p in range(len(patterns))
        for f in finish.keys()
    }
    m.param["a"] = a
    m.param["length_s"] = {s: stocks[s]["length"] for s in stocks.keys()}
    m.param["length_f"] = {f: finish[f]["length"] for f in finish.keys()}
    m.param["demand_finish"] = {f: finish[f]["demand"] for f in finish.keys()}

    ap_upper_bound = {
        f: max([int(stocks[s]["length"] / finish[f]["length"]) for s in stocks.keys()])
        for f in finish.keys()
    }
    m.param["ap_upper_bound"] = ap_upper_bound

    m.get_output("solve;")

    bp_values = dict(m.var["bp"].get_values())
    ap_values = dict(m.var["ap"].get_values())

    # Retrieve the patterns
    new_pattern = {
        "stock": [s for s in stocks.keys() if bp_values[s] > 0.5][0],
        "cuts": {f: round(ap_values[f]) for f in finish.keys()},
    }

    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': 'C', 'cuts': {'S': 1, 'M': 1, '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_{s_p}\) are relaxed to non-negative reals.

\[\begin{split} \begin{align} \min\quad & \sum_{p\in P} c_{s_p} x_{s_p} \\ \text{s.t.}\quad & \sum_{p\in P}a_{pf} x_{s_p} \geq d_f && \forall f\in F\\ & x_{s_p} \in \mathbb{R}_+ && \forall p\in P\\ \end{align} \end{split}\]

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\).

\[\begin{split} \begin{align} V_s = \max\quad & \left(\sum_{f\in F} \pi_{f} a'_{sf}\right) - c_s \\ \text{s.t.}\quad & \sum_{f\in F}l^F_f a'_{sf} \leq l^S_s && \\ & a'_{sf} \in \mathbb{Z}_+ && \forall f\in F\\ \end{align} \end{split}\]

The pattern demonstrating the largest return \(V_s\) is returned as a candidate to add the list of patterns.

This dual optimization problem might be implemented with two AMPL objects, each one representing a model. After declaring the model, it is only necessary to reset the data and assign it with new patterns or dual values. The “new_pattern_prob” corresponds to the problem with dual values, and “generate_pattern_dual_prob” to the relaxed one.

new_pattern_prob = AMPL()
new_pattern_prob.option["solver"] = SOLVER_MILO
new_pattern_prob.eval(
    """
        set F;

        param c;
        param ap_upper_bound{F};
        param demand_dual{F};

        # length stock
        param length_s;
        # length finished pieces
        param length_f{F};

        # Define second problems with new pattern in stock s
        var ap{f in F} integer >= 0 <= ap_upper_bound[f];

        maximize marginal_cost:
           sum{f in F} ap[f] * demand_dual[f] - c;

        subject to stock_length:
           sum{f in F} ap[f] * length_f[f] <= length_s;
"""
)


def new_pattern_problem(finish, length_s, cost_s, ap_upper_bound, demand_duals):
    m = new_pattern_prob
    m.eval("reset data;")
    m.set["F"] = list(finish.keys())
    m.param["c"] = cost_s
    m.param["length_s"] = length_s
    m.param["length_f"] = {f: finish[f]["length"] for f in finish.keys()}
    m.param["ap_upper_bound"] = ap_upper_bound
    m.param["demand_dual"] = demand_duals
    m.get_output("solve;")

    marg_cost = m.obj["marginal_cost"].value()
    pattern = {f: round(m.var["ap"][f].value()) for f in finish.keys()}
    return marg_cost, pattern


generate_pattern_dual_prob = AMPL()
generate_pattern_dual_prob.option["solver"] = SOLVER_MILO
generate_pattern_dual_prob.eval(
    """
        set F;
        set P;

        param c{P};
        param demand_finish{F};
        # how many F pieces are returned from pattern p
        param a{F, P};

        # Define first problem with known patterns
        var x{P} >= 0; # relaxed integrality

        minimize cost:
            sum{p in P} c[p] * x[p];

        subject to demand{f in F}:
           sum{p in P} a[f,p]*x[p] >= demand_finish[f];
"""
)


def generate_pattern_dual(stocks, finish, patterns):
    m = generate_pattern_dual_prob

    m.set["F"] = list(finish.keys())
    m.set["P"] = list(range(len(patterns)))

    s = {p: patterns[p]["stock"] for p in range(len(patterns))}
    c = {p: stocks[s[p]]["cost"] for p in range(len(patterns))}
    m.param["c"] = c

    a = {
        (f, p): patterns[p]["cuts"][f]
        for p in range(len(patterns))
        for f in finish.keys()
    }
    m.param["a"] = a
    m.param["demand_finish"] = {f: finish[f]["demand"] for f in finish.keys()}
    m.get_output("solve;")

    dual_values = dict(m.getConstraint("demand").get_values(suffixes="dual"))

    ap_upper_bound = {
        f: max([int(stocks[s]["length"] / finish[f]["length"]) for s in stocks.keys()])
        for f in finish.keys()
    }
    demand_duals = {f: dual_values[f] for f in finish.keys()}
    # use get_values() for getting the duals

    marginal_values = {}
    pattern = {}
    for s in stocks.keys():
        marginal_values[s], pattern[s] = new_pattern_problem(
            finish, stocks[s]["length"], stocks[s]["cost"], ap_upper_bound, demand_duals
        )

    s = max(marginal_values, key=marginal_values.get)
    new_pattern = {"stock": s, "cuts": pattern[s]}
    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_dual(stocks, finish, patterns)
{'stock': 'C', 'cuts': {'S': 1, 'M': 1, 'L': 1}}

The following cell compares the time required to generate a new pattern by the two methods for a somewhat larger data set.

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 generate_pattern_bilinear(stocks, finish, patterns)

print("Testing generate_patterns_dual: ", end="")
%timeit generate_pattern_dual(stocks, finish, patterns)
Testing generate_patterns_bilinear: 104 ms ± 36.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Testing generate_patterns_dual: 111 ms ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

A hybrid solution algorithm using pattern generation#

The following cell incorporates the two methods of pattern generation into a hybrid algorithm to solve the cutting stock problem. The algorithm starts with with make_patterns to create an initial list of patterns, then uses the linear dual to adds patterns until no new patterns are found. In phase two of the algorithm, the bilinear model is used to find additional patterns, if any, that may further reduce the objective function. There has been no exhaustive testing or attempt to compare this empirical method with others.

def cut_stock(stocks, finish):
    # Generate initial set of patterns
    patterns = make_patterns(stocks, finish)

    # Phase 1: Generate patterns using dual method
    print("Phase 1 ", end=".")
    new_pattern = generate_pattern_dual(stocks, finish, patterns)
    while new_pattern not in patterns:
        patterns.append(new_pattern)
        new_pattern = generate_pattern_dual(stocks, finish, patterns)
        print(end=".")

    x, cost = cut_patterns(stocks, finish, patterns)
    print(f" Cost = {cost}")

    # Phase 2: Generate patterns using bilinear method
    print("Phase 2 ", end=".")
    new_pattern = generate_pattern_bilinear(stocks, finish, patterns)
    while new_pattern not in patterns:
        patterns.append(new_pattern)
        new_pattern = generate_pattern_bilinear(stocks, finish, patterns)
        print(end=".")

    x, cost = cut_patterns(stocks, finish, patterns)
    print(f" Cost = {cost}")

    # Get the indices of non-zero patterns
    non_zero_indices = [index for index, value in enumerate(x) if value > 0]

    # Return only the non-zero patterns, their corresponding values, and the cost
    return (
        [patterns[index] for index in non_zero_indices],
        [x[index] for index in non_zero_indices],
        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.0
Phase 2 ... Cost = 170.0
<Axes: title={'center': 'Cost = 170.0'}>
../../_images/cb18734198a210b7f11dfbbb8e81f158faa36c37f6ed3241dc9154d41b7a23b0.png

Examples#

Example from JuMP documentation for column generation#

https://jump.dev/JuMP.jl/stable/tutorials/algorithms/cutting_stock_column_generation/#:~:text=The cutting stock problem is,while maximizing the total profit.

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 = 360.0
Phase 2 .. Cost = 360.0
<Axes: title={'center': 'Cost = 360.0'}>
../../_images/396e0f7cec191f69f0c16f0a65f6f7ae92a3c4d12f9c6852e900e5c7befe7427.png

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.0
Phase 2 .. Cost = 73.0
<Axes: title={'center': 'Cost = 73.0'}>
../../_images/663aba3f46a3fab384b3c584d3543a983cd06e56d4e596bab44f90858c422579.png

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

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")
# stocks = stocks.drop(["price"], axis=1)
if not "price" in stocks_df.columns:
    stocks["price"] = stocks_df["length"]
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 = tuple(set(finish_df["kind"]))

kerf = 0.25

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"] + 0 * kerf,
            "cost": stocks_df.loc[i, "price"],
        }

    patterns, x, cost = cut_stock(stocks, finish)
    plot_nonzero_patterns(stocks, finish, patterns, x, cost)
Kind = 2x4
Phase 1 ..... Cost = 32.440000000000005
Phase 2 .. Cost = 32.440000000000005
../../_images/e05fb713fd06ac90b90db34cc577e0323c22123541336637332add1c008f22c9.png

Purchase List

df = pd.DataFrame(patterns)
df["cuts"] = x
df["stock"] = df["stock"].astype(str)
df = df.pivot_table(index="stock", values="cuts", aggfunc="sum")
df.index = df.index.astype(int)
df = df.sort_index()

df
cuts
stock
84 2.0
144 5.0
168 1.0

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 Gamory 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