Course policies and submitting this assignment¶
Read the Artificial Intelligence Policy and Collaboration Policy and Honor Code before starting.
Submit two files to Canvas:
One scanned PDF containing all handwritten work for Problem 1 and 3-A. Put the work in problem order and label each part. A clear phone photo assembled into a single PDF is fine.
A copy of this notebook, with your code cells run.
The three AI categories are:
No AI.
AI permitted after independent work --- spend the stated time on your own first, then AI and genuine collaboration (including coding together) are permitted. Everyone must still contribute intellectually, understand the work, and verify it.
AI required --- independent work first, then using an AI tool is part of the problem.
Each problem below states its own category and independent-work window --- see the note at the start of each problem rather than a single rule for the whole assignment.
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 assignment has three problems:
formulate reactor-selection logic and reformulate it by hand;
compare big- and convex-hull formulations for strip packing; and
audit an AI-generated finite-difference model of transient heat conduction.
Problems 1 and 2 continue the thread that began in Problem 3 of Pyomo Homework 2. Problem 3 asks you to make an independent modeling plan before using AI, then verify the result with equations, structural diagnostics, and a mesh-refinement check.
| # | Problem | Format |
|---|---|---|
| 1 | Reactor-selection GDP | pencil and paper, no solver, no code cell |
| 2 | Strip packing: big- vs. convex hull | notebook/code |
| 3-A | Derive and plan the heat-conduction model | pencil and paper, no AI |
| 3-B -- 3-E | Implement, audit, and verify the model | notebook/code, AI required after 3-A |
# 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¶
AI category: AI permitted after independent work. Spend up to 30 minutes on paper, without AI and without the solution pages, before checking any help against it. 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 |
|---|---|
| low-conversion reactor (cheap) | |
| high-conversion reactor (expensive) | |
| absorber | |
| membrane | |
| cryogenic separation |
Write for the proposition “unit is installed”, and let 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 [million USD]:
| Reactor | Conversion window | [million USD] |
|---|---|---|
| 40 | ||
| 90 |
Conversion is physically bounded by .
Instruction. Same notation rules as Problem 3 of Pyomo Homework 2: capital for Booleans, lower-case for binaries, inequalities in the standard form .
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 and its binary .
1-B. Write each design rule as a logic proposition¶
Write rules 1 through 6 in the symbols using , , , and . 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 .
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 over to reach conjunctive normal form,
then substitute for and for , 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 .
Then answer: which of your six rules from 1-B is now redundant, because already says it?
1-E. Problem size¶
Let be the number of reactor options and the number of separator options. Using these symbols, determine the size of the model after the disjunction has been reformulated with big-. 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 and , 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- versus convex hull in Pyomo¶
AI category: AI permitted after independent work. Spend up to 30 minutes attempting 2-A on your own before using AI or the class notebook.
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 . Rectangle has length (along the strip) and height (across it), and is placed by the coordinates of its lower-left corner. The objective is to minimize the length of strip used, :
Non-overlap is a disjunction for every pair of rectangles with , with four disjuncts --- is left of , is right of , is below , or is above :
Nothing about the method changes from Problem 3 of Pyomo Homework 2 --- there are just disjunctions instead of one, and four terms instead of three.
Reference. Strip packing instance from the MINLP library,
https://
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- 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- 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 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 area bound from the note above. Then answer, in a few sentences each:
Do the two reformulations reach the same optimal ? 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 with width and height , plus a line marking the strip length . 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.bigmintroduced and say where the number comes from.The number of disjunctions grows as 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- 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.
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 3. Audit an AI-generated heat-conduction model¶
AI category: AI required after independent work. Spend 20 minutes deriving the model and planning its discretization without AI or solution code (3-A). Record that work on paper. Then ask an AI tool to help implement the model and audit its response (3-B onward).
Consider transient radial conduction in a slab of half-thickness :
with , symmetry at , and a prescribed surface temperature . Use dimensionless variables so that , , , and .
You may consult the course heat-conduction notebook and the ND Pyomo Cookbook chapter after recording your independent plan in 3-A. These are verification sources, so inspecting their code is allowed for this problem. This differs from Project 1, where the published implementation is intentionally withheld until after your independent validation.
3-A. Derive and plan without AI¶
Turn in a concise handwritten response. If any of these terms are unfamiliar, the linked pages are a fine place to start --- this problem is as much about learning to navigate PDE vocabulary you have not seen before as it is about the specific slab.
Write the initial condition and both boundary conditions.
Explain why the center condition is a Neumann condition and the surface condition is a Dirichlet condition.
Sketch a finite-difference stencil for the interior second derivative.
Predict what can go wrong if a central-difference transformation is applied before the center derivative is defined.
3-B. Generate a candidate implementation¶
Ask an AI tool to write a Pyomo.DAE implementation that:
discretizes space with finite differences and time with orthogonal collocation;
enforces the initial and boundary conditions without duplicate corner equations;
has zero degrees of freedom after discretization; and
solves with Ipopt.
Run the code, but do not treat a successful solve as verification. Save only the final code you audit; you do not need to submit a transcript.
3-C. Audit structure before solving¶
Ask an AI tool to help you read the output of this step --- the diagnostics below are dense, and
understanding what they mean is the point, not just running them. Two notebooks work through the same
kind of audit in more depth if you want another worked example: the course’s
heat-conduction notebook (this exact model, three geometries) and
NLP Diagnostics (build_model_size_report and DiagnosticsToolbox on
other models).
Use build_model_size_report() and DiagnosticsToolbox to check the generated model. Then answer:
Is the discretized model square after fixed variables are accounted for?
Which constraint defines the center derivative after spatial discretization?
Are the initial and surface conditions both imposed at the corner ? Explain how the implementation avoids a conflict.
Identify and correct at least one modeling or implementation weakness. If the first AI response is correct, perturb one boundary-condition line deliberately and show that your audit catches the defect.
# Add your solution here3-D. Verify the numerical result¶
Solve the corrected model on at least two spatial meshes. Report the center temperature at and the change caused by mesh refinement. Plot the temperature profile for the finer mesh.
Also compare the qualitative behavior with the course notebook:
temperature remains between the initial and surface temperatures;
the profile is symmetric about the slab center; and
heat propagates inward from the surface.
A visually plausible curve is evidence, but it is not a substitute for the structural audit.
# Add your solution here3-E. Compare against the analytic solution¶
3-D verified your model the way you would for a PDE with no known closed-form solution: mesh refinement and qualitative physical checks. This particular slab problem happens to have an exact series solution, so use it to see how good those proxies actually were.
Evaluate this series at the center (, ); a few dozen terms is plenty.
Report the absolute error between it and your center temperature from 3-D, for both meshes.
Does the error get smaller as you refine the mesh? Does it shrink smoothly, or level off?
Most PDEs you discretize in practice will not have a closed-form solution. Which of your 3-D checks would you still trust if this series did not exist?
# Add your solution here
Problem 3 AI use report¶
At the end of Problem 3, report for the problem as a whole:
approximately how long you spent on the 20-minute independent attempt in 3-A;
the result of that independent attempt, and how far you got;
what the AI tool produced or changed from 3-B onward;
the checks that established whether the result was trustworthy;
at least one defect you found or deliberately introduced in 3-C; and
one limitation of your verification.
If you used no AI, say so explicitly. Do not submit prompts or transcripts.