Extra material: Strip packing

Strip packing (SP) refers to the problem of packing rectangles onto a two dimensional strip of fixed width.

Many variants of this problem have been studied, the most basic is the pack a set of rectangles onto the shortest possible strip without rotation. Other variants allow rotation of the rectangles, require a packing that allows cutting the rectangles out of the strip with edge-to-edge cuts (guillotine packing), extends the strip to three dimensions, or extends the packing to non-rectangular shapes.

The extensive study of strip packing problems is motivated by their many industrial applications including

• placement of macro cells in semiconductor layouts,

• wood and textile cutting operations,

• laying out workstations in manufacturing facilities,

• allocating communications bandwidth between two endpoints,

• planning and scheduling \(CO_2\) utilization for enhanced oil recovery,

• scheduling allocations of a common resource.

Finding optimal solutions to strip packing problems is combinatorially difficult. Strip packing belongs to a class of problems called “NP-hard” for which known solution algorithms require effort that grows exponentially with problem size. For that reason much research on strip packing has been directed towards practical heuristic algorithms for finding good, though not optimal, for industrial applications.

Here we develop AMPL models to find optimal solutions to smaller but economically relevant problems. We use the problem of packing boxes onto shortest possible shelf of fixed width.

# install dependencies and select solver
%pip install -q amplpy matplotlib

SOLVER = "highs"

from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Problem Statment

Imagine a collection of \(N\) boxes that are to placed on shelf. The shelf depth is \(D\), and the dimensions of the boxes are \((w_i, d_i)\) for \(i=0, \ldots, N-1\). The boxes can be rotated, if needed, to fit on the shelf. How wide of a shelf is needed?

We will start by creating a function to generate a table of \(N\) boxes. For concreteness, we assume the dimensions are in millimeters.

import random
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle

# random seed
random.seed(1842)


# generate boxes
def generate_boxes(N, max_width=200, max_depth=200):
    boxes = pd.DataFrame()
    boxes["w"] = [random.randint(0.2 * max_width, max_width) for i in range(N)]
    boxes["d"] = [random.randint(0.2 * max_depth, max_depth) for i in range(N)]
    return boxes


N = 8
boxes = generate_boxes(8)
display(boxes)

# set shelf width as a multiple of the deepest box
D = 2 * boxes["d"].max()
print("Shelf Depth = ", D)

	w	d
0	82	103
1	73	48
2	171	53
3	73	99
4	167	85
5	151	172
6	54	130
7	126	94

Shelf Depth =  344

A lower and upper bounds on shelf width

A lower bound on the shelf width is established by from the area required to place all boxes on the shelf.

\[W_{lb} = \frac{1}{D}\sum_{i=0}^{N-1} w_i d_i\]

An upper bound is established by aligning the boxes along the front of the shelf without rotation. To set the stage for later calculations, the position of the rectangle on the shelf is defined by bounding box \((x_{i,1}, y_{i,1})\) and \((x_{i,2}, y_{i,2})\) extending from the lower left corner to the upper right corner.

\[\begin{split} \begin{align*} x_{i,2} & = x_{i,1} + w_i \\ y_{i,2} & = y_{i,1} + d_i \\ \end{align*} \end{split}\]

An additional binary variable \(r_i\) designates whether the rectangle has been rotated. The following cell performs these calculations to create and display a data frame showing the bounding boxes.

def pack_boxes_V0(boxes):
    soln = boxes.copy()
    soln["x1"] = soln["w"].cumsum() - soln["w"]
    soln["x2"] = soln["w"].cumsum()
    soln["y1"] = 0
    soln["y2"] = soln["d"]
    soln["r"] = 0
    return soln


def show_boxes(soln, D):
    """Show bounding boxes on a diagram of the shelf."""
    fig, ax = plt.subplots(1, 1, figsize=(14, 4))
    for i, x, y, w, h, r in zip(
        soln.index, soln["x1"], soln["y1"], soln["w"], soln["d"], soln["r"]
    ):
        c = "g"
        if r:
            h, w = w, h
            c = "r"
        ax.add_patch(Rectangle((x, y), w, h, edgecolor="k", facecolor=c, alpha=0.6))
        xc = x + w / 2
        yc = y + h / 2
        ax.annotate(
            i, (xc, yc), color="w", weight="bold", fontsize=12, ha="center", va="center"
        )

    W_lb = (soln["w"] * soln["d"]).sum() / D
    ax.set_xlim(0, 1.1 * soln["w"].sum())
    ax.set_ylim(0, D * 1.1)
    ax.axhline(D, label="shelf depth $D$", lw=0.8)
    ax.axvline(W_lb, label="lower bound $W_{lb}$", color="g", lw=0.8)
    ax.axvline(soln["x2"].max(), label="shelf width $W$", color="r", lw=0.8)
    ax.fill_between([0, ax.get_xlim()[1]], [D, D], color="b", alpha=0.1)
    ax.set_title(f"shelf width $W$ = {soln['x2'].max():.0f}")
    ax.set_xlabel("width")
    ax.set_ylabel("depth")
    ax.set_aspect("equal")
    ax.legend(loc="upper right")


soln = pack_boxes_V0(boxes)
display(soln)
show_boxes(soln, D)

	w	d	x1	x2	y1	y2	r
0	82	103	0	82	0	103	0
1	73	48	82	155	0	48	0
2	171	53	155	326	0	53	0
3	73	99	326	399	0	99	0
4	167	85	399	566	0	85	0
5	151	172	566	717	0	172	0
6	54	130	717	771	0	130	0
7	126	94	771	897	0	94	0

Modeling Strategy

At this point the reader may have some ideas on how to efficiently pack boxes on the shelf. For example, one might start by placing the larger boxes to left edge of the shelf, then rotating and placing the smaller boxes with a goal of minimized the occupied with of the shelf.

Modeling for optimization takes a different approach. The strategy is to describe constraints that must be satisfied for any solution to the problem, then let the solver find the a choice for the decision variables that minimizes width. The constraints include:

• The bounding boxes must fit within the boundaries of the shelf, and to the left of vertical line drawn at \(x = W\).

• The boxes can be rotated 90 degrees.

• The boxes must not overlap in either the \(x\) or \(y\) dimensions.

Version 1: An AMPL model to line up the boxes

For this first AMPL model we look to reproduce a lineup of the boxes on the shelf. In the case the problem is to minimize \(W\) where

\[\begin{split} \begin{align*} & \min W \\ \text{subject to:}\qquad\qquad \\ x_{i, 2} & = x_{i, 1} + w_i & \forall i\\ x_{i, 2} & \leq W & \forall i\\ x_{i, 1}, x_{i, 2} & \geq 0 & \forall i \\ \\ [x_{i, 2} \leq x_{j,1}] & \veebar [ x_{j, 2} \leq x_{i, 1}] & \forall i < j \\ \\ \end{align*} \end{split}\]

This first model does not consider rotation or placement of the boxes in the \(y\) dimension, so those decisions are not included.

The disjunctive constraints specify relationships between \(x_{i,1}\) and \(x_{i,2}\) to prevent overlapping positions of boxes on the shelf. The disjunctions require that either that box \(i\) is to the left of box \(j\) or that box \(j\) is the left of box \(i\). This is specified as an exclusive or disjunction because both conditions can be true at the same time. This disjunctive relationship must hold for every pair of boxes that are different from each other, but specifying \(i\) doesn’t overlap with \(j\) assures \(j\) doesn’t overlap \(i\). So it is only necessary to specify disjunctions for all pairs \(i, j\) where \(i < j\).

The corresponding AMPL model is a direct implementation of this model. One feature of the implementation is the use of a set PAIRS to identify the disjunctions. Defining this set simplifies coding for the corresponding disjunction.

%%writefile pack_boxes_V1.mod

set BOXES;
set PAIRS within {BOXES,BOXES};

param w{BOXES};

param W_ub;

var W >= 0, <= W_ub;
var x1{BOXES} >= 0, <= W_ub;
var x2{BOXES} >= 0, <= W_ub;

minimize minimize_width: W;

s.t. bounding_box {i in BOXES}:
    x2[i] == x1[i] + w[i];
s.t. width {i in BOXES}:
    x2[i] <= W;
s.t. no_overlap {(i, j) in PAIRS}:
    (x2[i] <= x1[j] or x2[j] <= x1[i])
    and not
    (x2[i] <= x1[j] and x2[j] <= x1[i]);

Overwriting pack_boxes_V1.mod

def pack_boxes_V1(boxes):
    # a simple upper bound on shelf width
    W_ub = boxes["w"].sum()

    BOXES = list(boxes.index)
    PAIRS = [(i, j) for i in BOXES for j in BOXES if i < j]

    m = AMPL()
    m.read("pack_boxes_V1.mod")

    m.set["BOXES"] = BOXES
    m.set["PAIRS"] = PAIRS

    m.param["w"] = boxes["w"]
    m.param["W_ub"] = int(W_ub)

    m.option["solver"] = SOLVER
    m.solve()

    soln = boxes.copy()
    soln["x1"] = m.var["x1"].to_dict()
    soln["x2"] = m.var["x2"].to_dict()
    soln["y1"] = 0
    soln["y2"] = soln["y1"] + soln["d"]
    soln["r"] = 0
    return soln


soln = pack_boxes_V1(boxes)
display(soln)
show_boxes(soln, D)

HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 896.9999999
344603 simplex iterations
73081 branching nodes
absmipgap=0.000150937, relmipgap=1.68269e-07

	w	d	x1	x2	y1	y2	r
0	82	103	0.0	82.0	0	103	0
1	73	48	770.0	843.0	0	48	0
2	171	53	432.0	603.0	0	53	0
3	73	99	82.0	155.0	0	99	0
4	167	85	603.0	770.0	0	85	0
5	151	172	281.0	432.0	0	172	0
6	54	130	843.0	897.0	0	130	0
7	126	94	155.0	281.0	0	94	0

Version 2: Rotating boxes

Rotating the boxes is an option for packing the boxes more tightly on the shelf. The boxes can be placed either in their original orientation or in a rotated orientation. This introduces a second exclusive or disjunction to the model that determines the orientation of the bounding box. A boolean indicator variable \(r_i\) tracks which boxes were rotated which is used in the show_boxes function to show which boxes have been rotated.

\[\begin{split} \begin{align*} & \min W \\ \text{subject to:}\qquad\qquad \\ x_{i, 2} & \leq W & \forall i\\ x_{i, 1}, x_{i, 2} & \geq 0 & \forall i \\ y_{i, 1} & = 0 & \forall i \\ \\ [x_{i, 2} \leq x_{j,1}] & \veebar [ x_{j, 2} \leq x_{i, 1}] & \forall i < j \\ \\ \begin{bmatrix} \neg r_i \\ x_{i,2} = x_{i,1} + w_i\\ y_{i,2} = y_{i,1} + d_i\\ \end{bmatrix} & \veebar \begin{bmatrix} r_i \\ x_{i,2} = x_{i,1} + d_i\\ y_{i,2} = y_{i,1} + w_i\\ \end{bmatrix} & \forall i < j \end{align*} \end{split}\]

For this version of the model the boxes will be lined up against the edge of the shelf with \(y_{i,1} = 0\). Decision variables are now included in the model for rotation \(r\) to the \(y\) dimension of the bounding boxes.

%%writefile pack_boxes_V2.mod

set BOXES;
set PAIRS within {BOXES,BOXES};

param w{BOXES};
param d{BOXES};

param W_ub;

var W >= 0, <= W_ub;
var x1{BOXES} >= 0, <= W_ub;
var x2{BOXES} >= 0, <= W_ub;
var y1{BOXES} >= 0, <= W_ub;
var y2{BOXES} >= 0, <= W_ub;
var r{BOXES} binary;

minimize minimize_width: W;

s.t. width {i in BOXES}:
    x2[i] <= W;
s.t. yloc {i in BOXES}:
    y1[i] == 0;
    
s.t. rotate {i in BOXES}:
    r[i] == 0 ==> 
        (x2[i] == x1[i] + w[i] and y2[i] == y1[i] + d[i])
    else
        (x2[i] == x1[i] + d[i] and y2[i] == y1[i] + w[i]);

s.t. no_overlap {(i, j) in PAIRS}:
    (x2[i] <= x1[j] or x2[j] <= x1[i])
    and not
    (x2[i] <= x1[j] and x2[j] <= x1[i])
    ;

Overwriting pack_boxes_V2.mod

def pack_boxes_V2(boxes):
    W_ub = boxes["w"].sum()

    BOXES = list(boxes.index)
    PAIRS = [(i, j) for i in BOXES for j in BOXES if i < j]

    m = AMPL()
    m.read("pack_boxes_V2.mod")

    m.set["BOXES"] = BOXES
    m.set["PAIRS"] = PAIRS

    m.param["w"] = boxes["w"]
    m.param["d"] = boxes["d"]
    m.param["W_ub"] = int(W_ub)

    m.option["solver"] = SOLVER
    m.solve()

    soln = boxes.copy()
    soln["x1"] = m.var["x1"].to_dict()
    soln["x2"] = m.var["x2"].to_dict()
    soln["y1"] = m.var["y1"].to_dict()
    soln["y2"] = m.var["y2"].to_dict()
    soln["r"] = m.var["r"].to_dict()

    return soln


soln = pack_boxes_V2(boxes)
display(soln)
show_boxes(soln, D)

HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 639.999502
7.40492e+06 simplex iterations
1.05588e+06 branching nodes
absmipgap=0.0608127, relmipgap=9.502e-05

	w	d	x1	x2	y1	y2	r
0	82	103	1.859995e+02	267.999502	0.0	103.0	0.0
1	73	48	5.300000e+01	101.000000	0.0	73.0	1.0
2	171	53	1.421085e-14	53.000000	0.0	171.0	1.0
3	73	99	5.669995e+02	639.999502	0.0	99.0	0.0
4	167	85	1.009995e+02	185.999502	0.0	167.0	1.0
5	151	172	3.619995e+02	512.999502	0.0	172.0	0.0
6	54	130	5.129995e+02	566.999502	0.0	130.0	0.0
7	126	94	2.679995e+02	361.999502	0.0	126.0	1.0

Version 3: Placing and Rotating boxes in two dimensions

Obviously the packages can be packed closer together by allowing boxes to be stacked deeper into the shelf. New constraints are needed to maintain the bounding boxes within the shelf depth, and to avoid overlaps in the \(y\) dimension.

\[\begin{split} \begin{align*} & \min W \\ \text{subject to:}\qquad\qquad \\ x_{i, 2} & \leq W & \forall i\\ y_{i, 2} & \leq D & \forall i\\ x_{i, 1}, x_{i, 2} & \geq 0 & \forall i \\ y_{i, 1}, y_{i, 2} & \geq 0 & \forall i \\ \\ \begin{bmatrix} x_{i, 2} \leq x_{j,1} \\ \end{bmatrix} \veebar \begin{bmatrix} x_{j, 2} \leq x_{i, 1} \\ \end{bmatrix} & \veebar \begin{bmatrix} y_{i, 2} \leq y_{j, 1} \\ \end{bmatrix} \veebar \begin{bmatrix} y_{j, 2} \leq y_{i, 1} \\ \end{bmatrix} & \forall i < j \\ \\ \begin{bmatrix} \neg r_i \\ x_{i,2} = x_{i,1} + w_i\\ y_{i,2} = y_{i,1} + d_i\\ \end{bmatrix} & \veebar \begin{bmatrix} r_i \\ x_{i,2} = x_{i,1} + d_i\\ y_{i,2} = y_{i,1} + w_i\\ \end{bmatrix} & \forall i < j \end{align*} \end{split}\]

%%writefile pack_boxes_V3.mod

set BOXES;
set PAIRS within {BOXES,BOXES};

param w{BOXES};
param d{BOXES};
param D;

param W_ub;

var W >= 0, <= W_ub;
var x1{BOXES} >= 0, <= W_ub;
var x2{BOXES} >= 0, <= W_ub;
var y1{BOXES} >= 0, <= W_ub;
var y2{BOXES} >= 0, <= W_ub;
var r{BOXES} binary;

minimize minimize_width: W;

s.t. width {i in BOXES}:
    x2[i] <= W;
s.t. height {i in BOXES}:
    y2[i] <= D;
    
s.t. rotate {i in BOXES}:
    r[i] == 0 ==> (x2[i] == x1[i] + w[i] and y2[i] == y1[i] + d[i])
    else
        (x2[i] == x1[i] + d[i] and y2[i] == y1[i] + w[i]);

s.t. no_overlap {(i, j) in PAIRS}:
    (x2[i] <= x1[j] or x2[j] <= x1[i] or y2[i] <= y1[j] or y2[j] <= y1[i])
    and not
    (x2[i] <= x1[j] and x2[j] <= x1[i] and y2[i] <= y1[j] and y2[j] <= y1[i])
    ;

Overwriting pack_boxes_V3.mod

def pack_boxes_V3(boxes, D):
    W_ub = boxes["w"].sum()

    BOXES = list(boxes.index)
    PAIRS = [(i, j) for i in BOXES for j in BOXES if i < j]

    m = AMPL()
    m.read("pack_boxes_V3.mod")

    m.set["BOXES"] = BOXES
    m.set["PAIRS"] = PAIRS

    m.param["w"] = boxes["w"]
    m.param["d"] = boxes["d"]
    m.param["W_ub"] = int(W_ub)
    m.param["D"] = int(D)

    m.option["solver"] = SOLVER
    m.solve()

    r = m.var["r"].to_dict()
    r = {i: int(round(r[i])) for i in r}

    soln = boxes.copy()
    soln["x1"] = m.var["x1"].to_dict()
    soln["x2"] = m.var["x2"].to_dict()
    soln["y1"] = m.var["y1"].to_dict()
    soln["y2"] = m.var["y2"].to_dict()
    soln["r"] = r
    return soln


soln = pack_boxes_V3(boxes, D)
display(soln)
show_boxes(soln, D)

HiGHS 1.5.1:HiGHS 1.5.1: optimal solution; objective 265.999999
2.86135e+06 simplex iterations
140133 branching nodes
absmipgap=8.4e-05, relmipgap=3.15789e-07

	w	d	x1	x2	y1	y2	r
0	82	103	0.000000	82.000000	151.0	254.0	0
1	73	48	217.999999	265.999999	126.0	199.0	1
2	171	53	94.999999	265.999999	205.0	258.0	0
3	73	99	0.000000	99.000000	271.0	344.0	1
4	167	85	99.000000	265.999999	258.0	343.0	0
5	151	172	0.000000	172.000000	0.0	151.0	1
6	54	130	82.000000	212.000000	151.0	205.0	1
7	126	94	171.999999	265.999999	0.0	126.0	1

Advanced Topic: Symmetry Breaking

One of the issues in combinatorial problem is the challenge of symmetries. A symmetry is a situation where a change in solution configuration leaves the objective unchanged. Strip packing problems are especially susceptible to symmetries.

Symmetries can significantly increase the effort needed to find and and verify an optimal solution. Trespalacios and Grossmann recently presented modification to the strip packing problem to reduce the number of symmetries. This is described in the following paper and implemented in the accompanying AMPL model.

Trespalacios, F., & Grossmann, I. E. (2017). Symmetry breaking for generalized disjunctive programming formulation of the strip packing problem. Annals of Operations Research, 258(2), 747-759.

%%writefile pack_boxes_V4.mod

set BOXES;
set PAIRS within {BOXES,BOXES};

param w{BOXES};
param d{BOXES};
param D;

param W_ub;

var W >= 0, <= W_ub;
var x1{BOXES} >= 0, <= W_ub;
var x2{BOXES} >= 0, <= W_ub;
var y1{BOXES} >= 0, <= W_ub;
var y2{BOXES} >= 0, <= W_ub;
var r{BOXES} binary;

minimize minimize_width: W;

s.t. width {i in BOXES}:
    x2[i] <= W;
s.t. height {i in BOXES}:
    y2[i] <= D;
    
s.t. rotate {i in BOXES}:
    r[i] == 0 ==> (x2[i] == x1[i] + w[i] and y2[i] == y1[i] + d[i])
    else
        (x2[i] == x1[i] + d[i] and y2[i] == y1[i] + w[i])
    ;

s.t. no_overlap {(i, j) in PAIRS}:
    x2[i] <= x1[j] or
    x2[j] <= x1[i] or
    (y2[i] <= y1[j] and x2[i] >= x1[j] + 1 and x2[j] >= x1[i] + 1) or
    (y2[j] <= y1[i] and x2[i] >= x1[j] + 1 and x2[j] >= x1[i] + 1);

Overwriting pack_boxes_V4.mod

def pack_boxes_V4(boxes, D):
    W_ub = boxes["w"].sum()

    BOXES = list(boxes.index)
    PAIRS = [(i, j) for i in BOXES for j in BOXES if i < j]

    m = AMPL()
    m.read("pack_boxes_V4.mod")

    m.set["BOXES"] = BOXES
    m.set["PAIRS"] = PAIRS

    m.param["w"] = boxes["w"]
    m.param["d"] = boxes["d"]
    m.param["W_ub"] = int(W_ub)
    m.param["D"] = int(D)

    m.option["solver"] = SOLVER
    m.solve()

    r = m.var["r"].to_dict()
    r = {i: int(round(r[i])) for i in r}

    soln = boxes.copy()
    soln["x1"] = m.var["x1"].to_dict()
    soln["x2"] = m.var["x2"].to_dict()
    soln["y1"] = m.var["y1"].to_dict()
    soln["y2"] = m.var["y2"].to_dict()
    soln["r"] = r

    return soln


soln = pack_boxes_V4(boxes, D)
display(soln)
show_boxes(soln, D)

HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 266
1.4118e+06 simplex iterations
74606 branching nodes
absmipgap=0.00274, relmipgap=1.03008e-05

	w	d	x1	x2	y1	y2	r
0	82	103	178.0	260.0	0.0	103.0	0
1	73	48	0.0	48.0	53.0	126.0	1
2	171	53	0.0	171.0	0.0	53.0	0
3	73	99	167.0	266.0	258.0	331.0	1
4	167	85	0.0	167.0	258.0	343.0	0
5	151	172	94.0	266.0	107.0	258.0	1
6	54	130	48.0	178.0	53.0	107.0	1
7	126	94	0.0	94.0	126.0	252.0	1