Pyomo Homework 3

Course policies and AI category¶

Read the Artificial Intelligence Policy and Collaboration Policy and Honor Code before starting. Assignment-specific directions control. The categories are No AI, AI permitted after independent work, and AI required.

Unless a problem says otherwise, its category is AI permitted after independent work. Spend about 30 minutes on each top-level problem without AI, solution pages, or another person’s help, stopping early if complete. You may consult lecture notes, textbooks, and nonsolution pages of the course website; bias toward those course sources. Afterward, AI and genuine collaboration, including coding together, are permitted. Everyone must contribute intellectually, understand the work, and verify it.

At the end of each top-level problem, add a concise AI and independent-work report: approximately how long the independent attempt took, how far you got, where you became stuck, any AI or collaborative help used afterward, and how you verified it. If you used no AI, say so. Do not submit prompts or transcripts. Time estimates help the instructor improve the assignment and are not a speed test.

How this assignment is organized¶

This homework has two problems, both from Logical Modeling and Generalized Disjunctive Programming. Problem 1 is pencil and paper --- no solver, no code cell --- and is written in the format of the in-person midterm. Problem 2 is code.

#	Problem	Format
1	Generalized disjunctive programming (GDP) modeling	pencil and paper
2	Strip packing: big- $M$ versus convex hull in Pyomo	code

Both problems continue the thread that began in Problem 3 of Pyomo Homework 2, where you reformulated a single disjunction by hand. Problem 1 asks you to write a disjunction from an English specification; Problem 2 hands both reformulations to Pyomo on a problem that is not small, and asks what each one costs. Several parts refer back to your answers in Pyomo Homework 2 --- keep it to hand.

Important

AI category: AI permitted after independent work. For each top-level problem, work for about 30 minutes without AI, solution pages, or help from another person, or stop early if you complete the problem. You may consult lecture notes, textbooks, and nonsolution pages of the course website; bias toward those course sources. Record your attempt before using AI. You may then use AI or collaborate with classmates, but you must understand and verify the work. Complete the short AI use report at the end of every problem.

In each report, state approximately how long the independent attempt took, how far you got, where you became stuck, any AI or collaborative help used afterward, and how you verified it. If you used no AI, say so. Use concise bullets; do not submit prompts or transcripts. Time estimates give the instructor useful data for improving the assignment and are not a speed test.

# Import the libraries you need here for the assignment

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

import pyomo.environ as pyo
import matplotlib.pyplot as plt

Problem 1. Generalized disjunctive programming (GDP) modeling¶

Pencil and paper. No solver, no code cell. This problem is written in the format of the in-person midterm.

Background. You are screening a superstructure for a new gas-processing plant. The design has two reactor options and three separator options:

Symbol	Unit
$R_1$	low-conversion reactor (cheap)
$R_2$	high-conversion reactor (expensive)
$A$	absorber
$M$	membrane
$CS$	cryogenic separation

Write $P_u$ for the proposition “unit $u$ is installed”, and let $y_u$ be the corresponding binary variable.

The process engineers hand you six design rules, in English:

Exactly one reactor is installed.
At least one separator is installed.
Cryogenic separation is only worth its capital cost when the high-conversion reactor is installed.
The absorber and the membrane cannot both be installed --- they compete for the same plot space.
If the high-conversion reactor is installed, then cryogenic separation or the membrane (or both) must be installed --- that reactor’s effluent carries light ends the absorber alone cannot remove.
If cryogenic separation is installed and the absorber is not, then the membrane must be installed.

The reactor choice also fixes the achievable conversion window and the reactor capital cost $c_R$ [million USD]:

Reactor	Conversion window	$c_R$ [million USD]
$R_1$	$0.30 \leq X \leq 0.60$	40
$R_2$	$0.80 \leq X \leq 0.95$	90

Conversion is physically bounded by $0 \leq X \leq 1$ .

Instruction. Same notation rules as Problem 3 of Pyomo Homework 2: capital $Y$ for Booleans, lower-case $y$ for binaries, inequalities in the standard form $g(x) \leq 0$ .

1-A. Translate the design into a mathematical model using set notation similar to our in-class examples¶

Next to each set, parameter, and variable, write a few-word description.

Sets.

Parameters (Data). Hint: state the units of each one.

Variables. Say for each one whether it is continuous, binary, or Boolean --- and be explicit about the relationship between a Boolean $Y_u$ and its binary $y_u$ .

1-B. Write each design rule as a logic proposition¶

Write rules 1 through 6 in the symbols $P_{R_1}, P_{R_2}, P_A, P_M, P_{CS}$ using $\neg$ , $\wedge$ , $\vee$ , $\veebar$ and $\Rightarrow$ . Do not convert to constraints yet.

Watch the parentheses. Rules 5 and 6 are the two where a careless reading gives the wrong proposition.

1-C. Convert the propositions to linear constraints¶

Convert each of rules 1 through 6 into linear constraints on the binaries $y_u$ .

For rules 1, 2, 3 and 4 you may quote the translation table from lecture --- give the row you used.

For rules 5 and 6, show the full three-step derivation:

① replace the implication ② apply De Morgan ③ distribute $\vee$ over $\wedge$ to reach conjunctive normal form,

then substitute $y_u$ for $P_u$ and $1 - y_u$ for $\neg P_u$ , and write one constraint per clause.

Finally, sanity-check your rule 6 constraint by substituting a case where the antecedent fires and a case where it does not.

1-D. Write the reactor choice as a disjunction¶

The conversion window and the reactor cost are not logic on binaries --- they switch a block of constraints on a continuous variable in or out. Write the reactor choice in GDP standard form: the disjunction over the reactor options, with the conversion bounds and the cost assignment inside each disjunct, and the appropriate $\Omega(Y)$ .

Then answer: which of your six rules from 1-B is now redundant, because $\Omega(Y)$ already says it?

1-E. Problem size¶

Let $N_r$ be the number of reactor options and $N_s$ the number of separator options. Using these symbols, determine the size of the model after the disjunction has been reformulated with big- $M$ . After each header below, give a number with a brief justification.

Number of continuous variables:

Number of integer/discrete variables:

Number of equality constraints:

Number of inequality constraints:

Then answer the interesting part: which of these counts actually grow with $N_r$ and $N_s$ , and which do not? Look carefully at the six logic constraints before you answer.

1-F. Classify the problem¶

Is the reformulated model an LP, QP, NLP, MILP or MINLP? Justify your classification by referring to the objective, the constraints and the variable domains.

Then: is the problem convex? Be careful --- answer separately for the model itself and for its relaxation, and connect your answer to what you found in Problem 3-D of Pyomo Homework 2.

Problem 1 AI use report¶

At the end of Problem 1, report for the problem as a whole:

approximately how long you spent on the single 30-minute independent attempt for Problem 1;
how far you got during that independent attempt;
where you got stuck, if anywhere;
any AI or collaborative help you used afterward; and
how you verified that help.

If you used no AI, say so explicitly. A few concise bullets are sufficient; do not submit prompts or transcripts. The time estimate gives the instructor useful data for improving the assignment and is not a speed test. Apply the 30-minute guideline once to the entire problem, not separately to each subpart.

Problem 2. Strip packing: big- $M$ versus convex hull in Pyomo¶

Problem 3 of Pyomo Homework 2 and Problem 1 above asked you to reformulate a disjunction by hand, on a problem small enough to see all of. This problem hands both reformulations to Pyomo on a problem that is not small, and asks what each one costs.

Background. Eight rectangles must be packed, without rotation and without overlap, into a strip of fixed width $W = 10$ . Rectangle $i \in N$ has length $L_i$ (along the strip) and height $H_i$ (across it), and is placed by the coordinates $(x_i, y_i)$ of its lower-left corner. The objective is to minimize the length of strip used, $lt$ :

\min_{x, y, lt} \; lt \qquad \text{s.t.} \qquad lt \geq x_i + L_i \quad \forall i \in N

(1)

Non-overlap is a disjunction for every pair of rectangles $(i,j)$ with $i < j$ , with four disjuncts --- $i$ is left of $j$ , $i$ is right of $j$ , $i$ is below $j$ , or $i$ is above $j$ :

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

(2)

Nothing about the method changes from Problem 3 of Pyomo Homework 2 --- there are just $\binom{8}{2} = 28$ disjunctions instead of one, and four terms instead of three.

Reference. Strip packing instance from the MINLP library, https://www.minlp.org/library/problem/index.php?i=121&lib=GDP (Vecchietti and Grossmann). The class notebook Modeling Disjunctions through the Strip Packing Problem works this model in full --- attempt 2-A yourself before opening it.

2-A. Build the GDP model¶

Complete create_model() below by adding the non-overlap disjunctions. Use the @model.Disjunction decorator (this is the notation you practised in Problem 1-J of Pyomo Homework 2) indexed over model.overlap_pairs, and return the four disjuncts as a Python list of expressions.

Do not write big- $M$ or convex hull constraints by hand. The whole point of writing the disjunction is that the transformation is somebody else’s job.

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


def create_model():
    """Build the strip packing problem as a generalized disjunctive program.

    Returns:
        model: Pyomo model, with the no-overlap disjunctions attached but NOT
            yet reformulated into a MILP.
    """
    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 (the 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}
    )

    # Extent of each rectangle along the strip (the x 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")

    # Upper bound on length: every rectangle laid end to end
    model.max_length = Param(
        initialize=sum(model.rect_length[i] for i in model.rectangles)
    )

    ## Variables
    model.x = Var(
        model.rectangles,
        bounds=(0, model.max_length),
        doc="Rectangle corner position along the strip",
    )

    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 position across the strip"
    )

    model.strip_length = Var(domain=NonNegativeReals, doc="Length of strip required")

    # The 28 unordered pairs of rectangles
    model.overlap_pairs = Set(
        initialize=model.rectangles * model.rectangles,
        dimen=2,
        filter=lambda b, i, j: i < j,
        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")

    ## Add the no-overlap disjunctions here!

    # Add your solution here

    return model


print(f"{len(create_model().overlap_pairs)} pairs of rectangles, one disjunction each")

2-B. Measure the size of a model¶

Before solving anything, write a function that reports the size of a Pyomo model: how many continuous variables, how many binary variables, and how many active constraints it has.

This is the same count you did by hand in Problem 3-E of Pyomo Homework 2, now done by the software. Use model.component_data_objects(Var, active=True) and the same for Constraint, and test each variable with .is_binary() and .is_continuous().

def model_size(model):
    """Count the continuous variables, binary variables, and active constraints.

    Argument:
        model: a Pyomo model, AFTER a GDP transformation has been applied

    Returns:
        dict with keys "continuous", "binary", "constraints"
    """
    # Add your solution here
    return {"continuous": n_cont, "binary": n_bin, "constraints": n_con}

2-C. Solve with the big- $M$ reformulation¶

Apply TransformationFactory("gdp.bigm") to a fresh model, solve it with HiGHS, and print the optimal strip length and the placement of every rectangle.

Hint: the solver returns coordinates a few floating-point units off a whole number (6.999999999999998, -0.0). Print with a format like %.4g --- that is numerical noise, not geometry.

def solve_variant(transformation, tee=False):
    """Build, transform, and solve the strip packing model.

    Argument:
        transformation: "gdp.bigm" or "gdp.hull"

    Returns:
        (model, size) -- the solved model and its size dict
    """
    model = create_model()
    # Add your solution here
    size = model_size(model)
    results = SolverFactory("appsi_highs").solve(model, tee=tee)
    assert check_optimal_termination(results), (
        f"Solve failed: status={results.solver.status}, "
        f"termination={results.solver.termination_condition}"
    )
    return model, size


bigm_model, bigm_size = solve_variant("gdp.bigm")

print("big-M reformulation")
print(f"  optimal strip length lt = {value(bigm_model.total_length):.4g}")
for i in bigm_model.rectangles:
    # +0.0 turns the solver's -0.0 into 0.0; .4g hides floating-point noise
    # such as 6.999999999999998, which is arithmetic, not geometry.
    xi = value(bigm_model.x[i]) + 0.0
    yi = value(bigm_model.y[i]) + 0.0
    print(f"  rectangle {i} at ({xi:.4g}, {yi:.4g})")
print(" ", bigm_size)

2-D. Solve with the convex hull reformulation, and compare¶

Do the same with gdp.hull. Then build the comparison.

Warning: A solve time depends on the machine, the solver version, and what else the computer is doing, so a timing you print here is not a number your classmate can reproduce.

Compare the two on things that are properties of the model:

Size: continuous variables, binary variables, constraints, from model_size.
Tightness of the relaxation: relax every binary to $0 \leq y \leq 1$ with TransformationFactory("core.relax_integer_vars") and solve the resulting LP. Its objective is the bound branch and bound starts from at the root node.

Report a table with both, plus the MILP optimum and the $109/10 = 10.9$ area bound from the note above. Then answer, in a few sentences each:

Do the two reformulations reach the same optimal $lt$ ? Do they reach it by the same placement? (2-E will show you.)
Which relaxation is tighter, and by how much? Compare each root bound to the true optimum and to the area bound.
The hull model is several times larger. Which kind of variable did it add --- and why does that make the trade worth taking? Connect this to your answer in Problem 3-E of Pyomo Homework 2.

def root_relaxation_bound(transformation):
    """Objective of the LP relaxation: the bound branch and bound starts from.

    Fully deterministic and machine-independent, unlike a solve time.
    """
    model = create_model()
    TransformationFactory(transformation).apply_to(model)
    # Add your solution here
    results = SolverFactory("appsi_highs").solve(model)
    assert check_optimal_termination(results)
    return value(model.total_length)


hull_model, hull_size = solve_variant("gdp.hull")

print("convex hull reformulation")
print(f"  optimal strip length lt = {value(hull_model.total_length):.4g}")
for i in hull_model.rectangles:
    # +0.0 turns the solver's -0.0 into 0.0; .4g hides floating-point noise
    # such as 6.999999999999998, which is arithmetic, not geometry.
    xi = value(hull_model.x[i]) + 0.0
    yi = value(hull_model.y[i]) + 0.0
    print(f"  rectangle {i} at ({xi:.4g}, {yi:.4g})")
print(" ", hull_size)

area_bound = sum(
    value(bigm_model.rect_length[i]) * value(bigm_model.rect_width[i])
    for i in bigm_model.rectangles
) / value(bigm_model.strip_width)

rows = [
    ("continuous variables", bigm_size["continuous"], hull_size["continuous"]),
    ("binary variables", bigm_size["binary"], hull_size["binary"]),
    ("constraints", bigm_size["constraints"], hull_size["constraints"]),
    ("LP relaxation bound at the root", root_relaxation_bound("gdp.bigm"),
     root_relaxation_bound("gdp.hull")),
    ("MILP optimum", value(bigm_model.total_length), value(hull_model.total_length)),
]

print(f"\n{'':35s} {'gdp.bigm':>10s} {'gdp.hull':>10s}")
for label, a, b in rows:
    print(f"{label:35s} {a:10.6g} {b:10.6g}")
print(f"\nArea lower bound (no packing can beat this): {area_bound:.4g}")
print("Wall-clock time is deliberately NOT reported: it is machine dependent.")

2-E. Visualize the packing¶

A table of coordinates is not a packing. Write a function that draws the solved model: one rectangle per item, drawn at $(x_i, y_i)$ with width $L_i$ and height $H_i$ , plus a line marking the strip length $lt$ . Use it on both solutions, one above the other.

import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle


def plot_packing(model, title, ax):
    """Draw a solved strip packing.

    Greyscale safety: every rectangle carries its own index printed in the
    middle AND its own hatch pattern, so no rectangle is identified by colour
    alone. The fills are light tints so that the black label and the black
    hatch stay legible in a photocopy.
    """
    # Light tints of the house Okabe-Ito palette. They are TINTS on purpose:
    # a saturated fill swallows the black index label and the black hatch, and
    # those two are what make the figure work in a photocopy.
    fills = ["#FFFFFF", "#FFE7B3", "#CDE6F5", "#CCEBDD"]
    hatches = ["", "///", "\\\\\\", "...", "xxx", "|||", "---", "+++"]

    W = value(model.strip_width)
    lt = value(model.strip_length)

    # Add your solution here

    ax.axvline(lt, color="black", linestyle="--", linewidth=1.5)
    ax.text(lt + 0.15, W * 0.94, f"$lt$ = {lt:g}", ha="left", va="top")
    ax.axhline(0, color="black", linewidth=1.0)
    ax.axhline(W, color="black", linewidth=1.0)
    ax.set_xlim(-0.4, max(lt, 1) + 3.5)
    ax.set_ylim(-0.4, W + 0.4)
    ax.set_aspect("equal")
    ax.set_xlabel("position along the strip, $x$")
    ax.set_ylabel("across, $y$")
    ax.set_title(title)
    return ax


fig, axes = plt.subplots(2, 1, figsize=(8, 7))
plot_packing(bigm_model, "gdp.bigm", axes[0])
plot_packing(hull_model, "gdp.hull", axes[1])
fig.tight_layout()
plt.show()

2-F. Discussion¶

Two or three sentences each.

Both reformulations describe the same set of feasible packings, yet 2-E almost certainly shows you two different pictures. Explain how both can be correct.
How do the disjunctions affect the degree of freedom analysis? Count the binaries that gdp.bigm introduced and say where the number comes from.
The number of disjunctions grows as $\binom{N}{2}$ in the number of rectangles. Contrast this with the logic constraints in 1-E, which did not grow at all. What is the difference between the two situations?
You wrote the model once, as a disjunction, and got two MILPs from it for free. Name one thing that would have gone wrong if you had typed the big- $M$ constraints by hand instead --- and refer to your answer in Problem 3-B part 3 of Pyomo Homework 2.

Problem 2 AI use report¶

At the end of Problem 2, report for the problem as a whole:

approximately how long you spent on the single 30-minute independent attempt for Problem 2;
how far you got during that independent attempt;
where you got stuck, if anywhere;
any AI or collaborative help you used afterward; and
how you verified that help.

Course policies and AI category¶

How this assignment is organized¶

Problem 1. Generalized disjunctive programming (GDP) modeling¶

1-A. Translate the design into a mathematical model using set notation similar to our in-class examples¶

1-B. Write each design rule as a logic proposition¶

1-C. Convert the propositions to linear constraints¶

1-D. Write the reactor choice as a disjunction¶

1-E. Problem size¶

1-F. Classify the problem¶

Problem 1 AI use report¶

Problem 2. Strip packing: big-MMM versus convex hull in Pyomo¶

2-A. Build the GDP model¶

2-B. Measure the size of a model¶

2-C. Solve with the big-MMM reformulation¶

2-D. Solve with the convex hull reformulation, and compare¶

2-E. Visualize the packing¶

2-F. Discussion¶

Problem 2 AI use report¶

Problem 2. Strip packing: big- $M$ versus convex hull in Pyomo¶

2-C. Solve with the big- $M$ reformulation¶