Modeling Disjunctions through the Strip Packing Problem

Prepared by: Prof. Alexander Dowling (adowling@nd.edu), Hailey Lynch (hlynch@nd.edu, 2023)

Introduction and Learning Objectives¶

This notebook illustrates a more complicated example of generalized disjunctive programs. Students will practice applying concepts of Logical Modeling and Generalized Disjunctive Programs to the Strip Packing Problem. Critical thinking discussion questions will be included to connect concepts from CBE 60499.

See the following notebook for a primer on Logical Modeling and Generalized Disjunctive Programs: Logical Modeling and Generalized Disjunctive Programs

Import Modules¶

# Imports
import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

milp_solver = "appsi_highs"

Pyomo.GDP: Strip Packing Problem¶

We will be looking at the Strip Packing Problem using Pyomo.GDP as an example for modeling disjunctions.

Aldo Vecchietti (1) and Ignacio Grossmann (2)
(1) INGAR – Instituto de Desarrollo y Diseño – CONICET – UTN, Avellaneda 3657 – Santa Fe ‐ Argentina
(2) Carnegie Mellon University, 5000 Forbes Av. ‐ Pittsburgh, PA ‐ USA

Additional Information	Link
Reference	https://www.minlp.org/library/problem/index.php?i=121&lib=GDP
Description	https://www.minlp.org/problems/ver/156/over/Strip-Packing-Overview.pdf
Results	https://www.minlp.org/problems/ver/156/results/Strip-Packing-Results.pdf
Problem Formulation	Sawaya and Grossmann (2005), A cutting plane method for solving linear generalized disjunctive programming problems

Problem Statement¶

The objective in this problem consists of minimizing the length of the strip $lt$ and by representing every rectangle by its coordinates in the $(x,y)$ space such that no overlap occurs between rectangles.

$\text{min} \ lt$

subject to...

Thus, every rectangle $i \in N$ has length $L_{i}$ , height $H_{i}$ , and coordinates $(x_{i}, y_{i})$ , where the point of reference corresponds to the upper left corner of every rectangle.

\ lt \geq x_{i} + L_{i} \ \ \forall i \in N

(1)

By constraining every pair of rectangles $(i,j)$ where $(i,j \in N, i<j)$ that no overlap occurs, we obtain a series of disjunctions with four disjuncts each, where each disjunct represents the position of rectangle $i$ in relation to rectangle $j$ .

\begin{bmatrix} Y_{ij}^{1} \\ x_{i} + L_{i} \leq x_{j} \end{bmatrix} \lor \begin{bmatrix} Y_{ij}^{2} \\ x_{j} + L_{j} \leq x_{i} \end{bmatrix} \lor \begin{bmatrix} Y_{ij}^{3} \\ y_{i} - H_{i} \geq y_{j} \end{bmatrix} \lor \begin{bmatrix} Y_{ij}^{4} \\ y_{j} - H_{j} \geq y_{i} \end{bmatrix}

(2)

Note that the y-coordinate of every rectangle is bounded from above by the fixed width of the strip $W$ , and that the upper bound $UB_{i}$ , which in a best case scenario would correspond to the optimal value of $lt$ , is obtained using a bottom-left rectangle-placing heuristic and serves as an upper bound for the x-coordinate of every rectangle.

x_{i} \leq UB_{i} - L_{i} \ \ \forall i \in N \tag{4}

(3)

H_{i} \leq y_{i} \leq W \ \ \forall i \in N \tag{5}

(4)

lt, x_{i}, y_{i} \in \mathbb{R}_{+}^{1}, Y_{ij}^{1}, Y_{ij}^{2}, Y_{ij}^{3}, Y_{ij}^{4} \in \{\text{True, False}\} \ \ \forall i, j \in N, i < j

(5)

Define model in Pyomo with GDP¶

First we will define the model for the Strip Packing Problem in Pyomo.

"""
Instead of using
# import pyomo.environ as pyo
We can import specific functions/objects
"""

from pyomo.environ import (
    check_optimal_termination,
    ConcreteModel,
    NonNegativeReals,
    Objective,
    Param,
    Set,
    SolverFactory,
    TransformationFactory,
    Var,
    value,
)

# Strip-packing example from http://minlp.org/library/lib.php?lib=GDP

# This model packs a set of rectangles without rotation or overlap within a
# strip of a given width, minimizing the length of the strip.


def create_model():
    """
    Build the strip packing problem model.

    Return:
    model: Pyomo model

    """

    ## Model
    model = ConcreteModel(name="Rectangles strip packing")

    ## Sets
    model.rectangles = Set(ordered=True, initialize=[0, 1, 2, 3, 4, 5, 6, 7])

    ## Parameters
    # Extent of each rectangle across the width of the strip (y direction)
    model.rect_width = Param(
        model.rectangles, initialize={0: 3, 1: 3, 2: 2, 3: 2, 4: 3, 5: 5, 6: 7, 7: 7}
    )

    # Parameter indexed by each rectangle.
    # rect_length is the extent along the strip (the x direction);
    # rect_width is the extent across the strip (the y direction).
    model.rect_length = Param(
        model.rectangles, initialize={0: 4, 1: 3, 2: 2, 3: 2, 4: 3, 5: 3, 6: 4, 7: 4}
    )

    model.strip_width = Param(initialize=10, doc="Available width of the strip")

    # Upperbound on length (default is sum of lengths of rectangles)
    model.max_length = Param(
        initialize=sum(model.rect_length[i] for i in model.rectangles),
        doc="Maximum length of the strip (if all rectangles were arranged "
        "lengthwise)",
    )

    ## Variables
    # x (length) and y (width) coordinates of each of the rectangles
    model.x = Var(
        model.rectangles,
        bounds=(0, model.max_length),
        doc="Rectangle corner x-position (position across length)",
    )

    # Width bounds
    def w_bounds(b, i):
        return (0, b.strip_width - b.rect_width[i])

    model.y = Var(
        model.rectangles,
        bounds=w_bounds,
        doc="Rectangle corner y-position (position down width)",
    )
    # Strip length
    model.strip_length = Var(domain=NonNegativeReals, doc="Length of strip required.")

    # Rectangle conflicts
    def rec_pairs_filter(b, i, j):
        return i < j

    model.overlap_pairs = Set(
        initialize=model.rectangles * model.rectangles,
        dimen=2,
        filter=rec_pairs_filter,
        doc="Set of possible rectangle conflicts",
    )

    ## Constraints
    @model.Constraint(model.rectangles)
    def strip_ends_after_last_rec(b, i):
        return b.strip_length >= b.x[i] + b.rect_length[i]

    ## Objective
    model.total_length = Objective(expr=model.strip_length, doc="Minimize length")

    ## Insert the no-overlap disjunctions here!

    # Add your solution here

    return model

@model.Disjunction(
    model.overlap_pairs,
    doc="Make sure that none of the rectangles on the strip overlap in "
    "either the x or y dimensions.",
)
def no_overlap(b, i, j):
    return [
        b.x[i] + b.rect_length[i] <= b.x[j],
        b.x[j] + b.rect_length[j] <= b.x[i],
        b.y[i] + b.rect_width[i] <= b.y[j],
        b.y[j] + b.rect_width[j] <= b.y[i],
    ]

Transform and Solve with Big M Relaxation¶

Now we will create the model and use Big-M Relaxation.

Big-M Implementation in Pyomo¶

# Creating the model
model = create_model()

Next, let’s transform the model. The updated model is really big, so we will not print it. This is because the transformation factory replaced the disjunctions with many Big-M constraints.

# Applying Big-M relaxation to the model
# Add your solution here

Finally, we’ll solve the model and examine the solution.

# Solve and print the solution
results = SolverFactory(milp_solver).solve(model, tee=True)
assert check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)
for i in model.rectangles:
    # %.4g, not %s: the solver returns coordinates a few ULPs off a whole
    # number (e.g. 6.999999999999998, -0.0), which is noise, not geometry.
    xi, yi = value(model.x[i]) + 0.0, value(model.y[i]) + 0.0
    print(f"Rectangle {i}: ({xi:.4g}, {yi:.4g})")
model.total_length.display()

Rectangle 0: (7, 0)
Rectangle 1: (4, 0)
Rectangle 2: (0, 0)
Rectangle 3: (2, 0)
Rectangle 4: (0, 7)
Rectangle 5: (0, 2)
Rectangle 6: (7, 3)
Rectangle 7: (3, 3)
total_length : Size=1, Index=None, Active=True
    Key  : Active : Value
    None :   True :  11.0

Transform and Solve with Convex Hull Relaxation¶

We will repeat the procedure above using Convex Hull Relaxation.

Convex Hull Implementation in Pyomo¶

# Creating the model
model = create_model()

# Applying convex hull relaxation to the model
# Add your solution here

# Solve and print the solution
results = SolverFactory(milp_solver).solve(model, tee=True)
assert check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)
for i in model.rectangles:
    # %.4g, not %s: the solver returns coordinates a few ULPs off a whole
    # number (e.g. 6.999999999999998, -0.0), which is noise, not geometry.
    xi, yi = value(model.x[i]) + 0.0, value(model.y[i]) + 0.0
    print(f"Rectangle {i}: ({xi:.4g}, {yi:.4g})")
model.total_length.display()

Discussion Questions¶

What is the advantage of using the decorator notation for optimization problems?
How do disjunctions affect the degree of freedom analysis? For eight rectangles, how many binary variables does gdp.bigm create, and where does that number come from?
Why is it necessary to make separate variables for length and width?
Compare the outputs of Big-M and Convex Hull. How are they different (if at all)?