Algorithms Homework 6 - Optimization for Decision Science

Assignment overview. Notebook MINLP Algorithms implemented branch and bound, outer approximation and generalized Benders decomposition by hand, on problems small enough to watch. This assignment does the opposite: you drive a production MINLP solver — Pyomo’s MindtPy — on a problem too big to watch, and use the theory from the handout to explain what comes back.

You will not write an algorithm here. You will write one line, over and over, with one option changed each time:

results = pyo.SolverFactory("mindtpy").solve(
    m, strategy="OA", mip_solver="appsi_highs", nlp_solver="ipopt"
)

and then account for the iteration counts, the run times, and — in two cases — the wrong answers.

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.

Tips and tricks¶

MindtPy’s strategy= argument selects the algorithm. "OA" is outer approximation (handout §4), "ECP" is the extended cutting plane method, "GOA" is global outer approximation, and "FP" is the feasibility pump.
MindtPy is a meta-solver. It has no numerics of its own: it calls an MILP solver for the master problems (mip_solver=) and an NLP solver for the subproblems (nlp_solver=). Both must be named.
Check the termination condition every single time. This is not boilerplate in this assignment — it is one of the graded findings. At least one strategy below returns an answer and a termination condition that is not optimal, and a student who skips the check will report a number that is 10 % off as if it were the optimum.
results.solver.iterations holds the number of major iterations for most strategies. Use getattr(results.solver, "iterations", None) so a missing field does not crash your loop.
MindtPy is chatty. Wrap the solve in the quiet() helper provided below when you want a clean table.
Wall time is measured with time.perf_counter(), and it is noisy at this problem size. Report it to two decimals, do not over-interpret differences under about 0.5 s, and re-run before you claim anything about speed.

Reference: the MindtPy documentation, and Biegler, Grossmann & Westerberg (1997) Appendix A, §A.3.2 and §A.3.4.

# This code cell installs packages on Colab

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 contextlib
import io
import time

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyomo.environ as pyo

rng = np.random.default_rng(seed=0)


@contextlib.contextmanager
def quiet():
    """Swallow a solver's chatter so a loop can print a clean table."""
    buf = io.StringIO()
    with contextlib.redirect_stdout(buf):
        yield buf

Part 0. The model¶

A whole-building insulation retrofit. Four wall assemblies — three walls and a roof — are to be insulated from a catalogue of eight materials. This is the notebook’s single-wall problem with three things added, each of which makes the binary structure genuinely hard:

Per-assembly layer choice. $y_{w,n} = 1$ if material $n$ is installed on assembly $w$ , with thickness $x_{w,n} \le t^{\max}_n y_{w,n}$ .
A stocking decision. $z_n = 1$ if material $n$ is used anywhere. Each stocked material costs a fixed crew/logistics charge, so the optimizer is pushed toward using few materials on many walls — the constraint $y_{w,n} \le z_n$ couples the assemblies together.
A capital budget. Total installed cost is capped, which is what makes this a constrained combinatorial problem rather than four independent ones.

\begin{aligned} \min_{x,y,z} \quad & \sum_{w} A_w \left[ \frac{\alpha}{R_w} + \beta \sum_n \left( a_n y_{w,n} + b_n x_{w,n} \right) \right] + \beta \, \kappa \sum_n z_n \\ \text{s.t.} \quad & R_w = R_{0,w} + \sum_n \frac{x_{w,n}}{k_n} && \forall w \\ & x_{w,n} \le t^{\max}_n \, y_{w,n}, \qquad y_{w,n} \le z_n && \forall w, n \\ & \sum_n x_{w,n} \le t^{\max}_w && \forall w \\ & \sum_{w} A_w \sum_n \left( a_n y_{w,n} + b_n x_{w,n} \right) + \kappa \sum_n z_n \le \mathcal{B} \\ & x \ge 0, \qquad y, z \in \{0,1\} \end{aligned}

(1)

The objective is convex for fixed $y$ : $\alpha/R_w$ is convex in $x$ because $R_w$ is affine and positive, and everything else is linear. Every constraint is linear. So this is a convex MINLP, and outer approximation and generalized Benders both apply with their guarantees intact. Part 4 breaks that on purpose.

There are $4 \times 8 + 8 = 40$ binary variables here, so complete enumeration would need $2^{40} \approx 10^{12}$ NLP solves.

# k: thermal conductivity [W/m/K]; a: fixed installation cost [$/m^2];
# b: installed material cost [$/m^3]; tmax: largest available thickness [m]
materials = pd.DataFrame(
    {
        "Fiberglass batt": {"k": 0.040, "a": 4.0, "b": 60.0, "tmax": 0.06},
        "Mineral wool": {"k": 0.030, "a": 5.0, "b": 150.0, "tmax": 0.06},
        "Rigid foam (low R)": {"k": 0.030, "a": 8.0, "b": 120.0, "tmax": 0.05},
        "Rigid foam (high R)": {"k": 0.015, "a": 8.0, "b": 180.0, "tmax": 0.05},
        "Aerogel blanket": {"k": 0.013, "a": 12.0, "b": 900.0, "tmax": 0.02},
        "Cellulose fill": {"k": 0.038, "a": 3.0, "b": 45.0, "tmax": 0.08},
        "Polyiso board": {"k": 0.022, "a": 7.0, "b": 160.0, "tmax": 0.05},
        "Vacuum panel": {"k": 0.007, "a": 25.0, "b": 2500.0, "tmax": 0.015},
    }
).T

# area [m^2]; R0: resistance of the structural elements [m^2 K / W];
# tmax: cavity depth available on that assembly [m]
assemblies = pd.DataFrame(
    {
        "North wall": {"area": 40.0, "R0": 2.0, "tmax": 0.15},
        "East wall": {"area": 25.0, "R0": 1.6, "tmax": 0.12},
        "South wall": {"area": 40.0, "R0": 2.4, "tmax": 0.10},
        "Roof": {"area": 60.0, "R0": 3.0, "tmax": 0.25},
    }
).T

ALPHA = 60.0  # annualized energy cost per unit U   [$ K / W / m^2]
BETA = 0.05  # equivalent annual cost factor on capital
CREW = 150.0  # fixed crew/logistics charge per material stocked [$]
BUDGET = 900.0  # capital budget [$]

display(materials)
display(assemblies)

def build_retrofit_minlp(budget=BUDGET):
    """Build the whole-building insulation retrofit MINLP.

    Arguments:
        budget: capital budget on installed cost [$]

    Returns:
        a Pyomo ConcreteModel
    """
    m = pyo.ConcreteModel("insulation retrofit MINLP")
    m.N = pyo.Set(initialize=list(materials.index))
    m.W = pyo.Set(initialize=list(assemblies.index))

    m.x = pyo.Var(m.W, m.N, domain=pyo.NonNegativeReals, bounds=(0, 0.3))
    m.y = pyo.Var(m.W, m.N, domain=pyo.Binary)  # material n installed on assembly w
    m.z = pyo.Var(m.N, domain=pyo.Binary)  # material n stocked at all

    @m.Constraint(m.W, m.N)
    def layer_available(m, w, n):
        return m.x[w, n] <= materials.loc[n, "tmax"] * m.y[w, n]

    @m.Constraint(m.W, m.N)
    def stocked(m, w, n):
        return m.y[w, n] <= m.z[n]

    @m.Constraint(m.W)
    def thickness_budget(m, w):
        # [m] cavity depth available on that assembly
        return sum(m.x[w, n] for n in m.N) <= assemblies.loc[w, "tmax"]

    @m.Constraint()
    def capital_budget(m):
        # [$] installed cost, before annualization
        return (
            sum(
                assemblies.loc[w, "area"]
                * (
                    materials.loc[n, "a"] * m.y[w, n]
                    + materials.loc[n, "b"] * m.x[w, n]
                )
                for w in m.W
                for n in m.N
            )
            + CREW * sum(m.z[n] for n in m.N)
            <= budget
        )

    @m.Objective(sense=pyo.minimize)
    def cost(m):
        total = 0
        for w in m.W:
            # [m^2 K / W] series resistance of the assembly
            R = assemblies.loc[w, "R0"] + sum(
                m.x[w, n] / materials.loc[n, "k"] for n in m.N
            )
            total += assemblies.loc[w, "area"] * (
                ALPHA / R
                + BETA
                * sum(
                    materials.loc[n, "a"] * m.y[w, n]
                    + materials.loc[n, "b"] * m.x[w, n]
                    for n in m.N
                )
            )
        total += BETA * CREW * sum(m.z[n] for n in m.N)
        return total

    return m


def report(m, label):
    """Print the stocked materials and the layer schedule of a solved model."""
    print(f"{label}: annualized cost = {pyo.value(m.cost):.4f} $/yr")
    print(f"  stocked: {[n for n in m.N if pyo.value(m.z[n]) > 0.5]}")
    for w in m.W:
        layers = {
            n: round(pyo.value(m.x[w, n]), 4) for n in m.N if pyo.value(m.y[w, n]) > 0.5
        }
        print(f"  {w:12s} {layers}")

The reference answer¶

Before benchmarking anything, get an answer you trust. bonmin is a dedicated MINLP solver (it implements NLP-based branch and bound and outer approximation, among others) and it is exact for convex problems, which this one is. Every strategy in Parts 1–3 is measured against this number.

m_ref = build_retrofit_minlp()
t0 = time.perf_counter()
results = pyo.SolverFactory("bonmin").solve(m_ref)
t_ref = time.perf_counter() - t0

assert pyo.check_optimal_termination(
    results
), f"bonmin failed: {results.solver.termination_condition}"

Z_REF = pyo.value(m_ref.cost)
report(m_ref, "bonmin")
print(f"  wall time = {t_ref:.2f} s")

bonmin: annualized cost = 2764.3022 $/yr
  stocked: ['Cellulose fill']
  North wall   {'Cellulose fill': 0.08}
  East wall    {'Cellulose fill': 0.08}
  South wall   {'Cellulose fill': 0.08}
  Roof         {}
  wall time = 0.63 s

Problem 1. One solve, and one check (10 points)¶

Solve the same model with MindtPy’s outer approximation strategy, using appsi_highs for the master MILPs and ipopt for the NLP subproblems.

Then, before printing anything:

Assert that the termination condition is optimal.
Assert that the objective agrees with Z_REF to within 10^-3.

Print the objective, the number of major iterations, and the wall time.

m1 = build_retrofit_minlp()

# Add your solution here

Question 1a. Outer approximation solves one NLP and one MILP per major iteration. Look at the iteration count you just printed. How many NLP solves is that, and how does it compare with the 2⁴⁰ that complete enumeration would need? What is doing the work of eliminating the rest?

# Add your solution here

Problem 2. The strategy sweep (30 points)¶

MindtPy 6.10.1 advertises four strategies. Try all four on the same model and tabulate what happens.

Do not assume they all work. One of them requires a compiled library that is not part of this course environment. Your loop must catch the exception, record it, and keep going — a benchmark that dies on the third row is not a benchmark.

Write a function benchmark(strategy, **options) that

builds a fresh model (never reuse a solved one — the warm start would contaminate the comparison),
times the solve,
returns a dict with keys strategy, termination, objective, iterations, time [s], and gap vs bonmin,
and on an exception returns the exception type and message in place of the numbers.

Then run it for "OA", "ECP", "GOA", "FP" and display the table.

def benchmark(strategy, model_builder=build_retrofit_minlp, reference=None, **options):
    """Run one MindtPy strategy on a fresh model and summarize the outcome.

    Arguments:
        strategy: MindtPy strategy string, e.g. "OA"
        model_builder: zero-argument callable returning a fresh ConcreteModel
        reference: objective value to compare against, or None
        **options: extra keyword arguments passed to MindtPy

    Returns:
        a dict of results; on failure, the exception type and message
    """
    # Add your solution here


# Add your solution here

Question 2a — GOA. One row failed. Report the exception, then explain in one sentence what MindtPy was trying to do when it failed, and why global outer approximation needs machinery that outer approximation does not. The following cell will help you diagnose it.

from pyomo.contrib.mcpp.pyomo_mcpp import mcpp_available

print("MC++ library available:", mcpp_available())

MC++ library available: False

# Add your solution here

Question 2b — FP. The feasibility pump row is the trap this assignment is built around. Report its termination condition and its objective. Is the objective feasible? Is it optimal? Would pyo.check_optimal_termination have caught it?

Then explain, in terms of the two bounds $z^u$ and $z^L$ from the handout, why the feasibility pump can never return a proof of optimality no matter how long you run it.

# Add your solution here

Question 2c — ECP. The extended cutting plane method reaches an objective close to, but not exactly, bonmin’s, and it does so much faster than OA. Report both numbers. Then explain the mechanism: what does ECP do at each iteration that OA does not do, and how does that account for both the speed and the residual error?

# Add your solution here

Problem 3. Toggling outer approximation (25 points)¶

strategy="OA" is not one algorithm; it is a family. Three options change it substantially.

(a) init_strategy= — where the first $y^{1}$ comes from. In the notebook we started from a deliberately poor guess (aerogel alone) and watched it cost an iteration. MindtPy offers

value	what it does
`"rNLP"`	solve the fully relaxed NLP (integrality dropped) and round — the default
`"max_binary"`	maximize the sum of the binaries subject to the linear constraints
`"initial_binary"`	use the values already stored in the model’s binary variables
`"FP"`	run the feasibility pump from Problem 2 as an initializer

Benchmark all four. "initial_binary" requires that every binary already have a value — set them yourself first, and use the answer to Problem 1 as a warm start to see the best case.

(b) add_no_good_cuts= — whether to add the integer cut that forbids revisiting a $y^{k}$ . Try True and False.

(c) single_tree= — MindtPy’s LP/NLP-based branch and bound, which solves one MILP tree and calls the NLP from a solver callback instead of restarting the tree every iteration. Try it and report what happens.

# Add your solution here

Question 3a. Every row that ran returned the same objective. Given that, what is the initialization buying, and which row bought the most? Tie your answer to the role $y^{1}$ plays in the handout’s flowchart.

# Add your solution here

Question 3b. single_tree=True failed. Report the message, and explain why the LP/NLP-based branch and bound method cannot be implemented on top of appsi_highs. What would you have to install?

# Add your solution here

Question 3c. add_no_good_cuts=False still converged, and to the right answer. In the notebook’s hand-written outer_approximation, removing the integer cut would have caused an infinite loop. Why does MindtPy survive without it?

# Add your solution here

Problem 4. Where outer approximation stops being an algorithm (25 points)¶

Everything above rested on convexity. The handout is explicit about it: the tangent plane in the master problem is a valid underestimator only because $f$ and $g$ are convex, and without that the “lower bound” is not a bound at all.

Here is the notebook’s five-material single-wall problem with one change: the capital cost of a layer has economies of scale,

\text{capital}_n(x_n) \;=\; b_n \, s \left[ (x_n + \epsilon)^{p} - \epsilon^{p} \right], \qquad p = 0.6, \quad s = 8, \quad \epsilon = 10^{-4},

(2)

which is concave in $x_n$ — the second unit of thickness costs less than the first, because the crew is already on site. That is entirely realistic, and it makes the objective nonconvex.

SCALE = 8.0  # cost multiplier on the concave capital term
EXPO = 0.6  # < 1 makes it concave: economies of scale
EPS = 1e-4  # offset keeping the derivative finite at x = 0
R0 = 2.0  # resistance of the structural elements  [m^2 K / W]
T_TOTAL = 0.15  # total thickness the wall cavity allows [m]
NC_MATERIALS = list(materials.index[:5])


def build_nonconvex_minlp(x_init=0.005):
    """Single-wall insulation MINLP with a concave (economies-of-scale) capital cost.

    Arguments:
        x_init: value every thickness variable is initialized at [m]

    Returns:
        a Pyomo ConcreteModel -- NONCONVEX
    """
    m = pyo.ConcreteModel("nonconvex insulation MINLP")
    m.N = pyo.Set(initialize=NC_MATERIALS)
    m.x = pyo.Var(m.N, domain=pyo.NonNegativeReals, bounds=(0, 0.06), initialize=x_init)
    m.y = pyo.Var(m.N, domain=pyo.Binary, initialize=0)

    @m.Constraint(m.N)
    def layer_available(m, n):
        return m.x[n] <= materials.loc[n, "tmax"] * m.y[n]

    @m.Constraint()
    def thickness_budget(m):
        return sum(m.x[n] for n in m.N) <= T_TOTAL

    @m.Objective(sense=pyo.minimize)
    def cost(m):
        R = R0 + sum(m.x[n] / materials.loc[n, "k"] for n in m.N)
        return ALPHA / R + BETA * sum(
            materials.loc[n, "a"] * m.y[n]
            + materials.loc[n, "b"] * SCALE * ((m.x[n] + EPS) ** EXPO - EPS**EXPO)
            for n in m.N
        )

    return m

(a) Plot the capital cost of one material against thickness, both the original linear form $b_n x_n$ and the concave form, and confirm visually that the new one is concave.

(b) Solve the model with couenne, which is a spatial branch and bound code and is therefore a rigorous global solver — it does not assume convexity. This is the reference.

(c) Solve it with MindtPy strategy="OA" from six different initializations (x_init $\in \{0, 0.001, 0.005, 0.01, 0.03, 0.06\}$ ) and tabulate the objective, the design, and the termination condition for each.

# Add your solution here

# Add your solution here

# Add your solution here

Question 4a. Summarize the table in two sentences. What is the largest error, in percent? What termination condition did MindtPy report in the rows that were wrong?

# Add your solution here

Question 4b. Explain the failure mechanism precisely, in the language of the handout. Which of the two bounds is invalid, and what does the algorithm’s termination test then actually test? Draw the tangent plane to the concave cost at $x = 0.005$ on your plot from (a) if it helps.

# Add your solution here

Question 4c. What would be the right tool here, and what would it cost? Name three options, and for each say what it gives up. (One of them you already met in Problem 2.)

# Add your solution here

Problem 5. Short answer (10 points)¶

Answer in two or three sentences each. These are the connections the assignment is testing; no code is needed.

5a. The notebook’s generalized Benders implementation took four major iterations where outer approximation took three, on the same problem from the same start. MindtPy has no strategy="GBD". Given what you found in Problem 2 about the cost of each strategy, and the bound inequality $z^L_{\text{GBD}} \le z^L_{\text{OA}}$ , why might that omission be defensible — and when would you want it back?

# Add your solution here

5b. Handout §2 says pure cutting-plane methods are not used to solve integer programs because of numerical instability, and §3 says the fix is branch and cut. ECP in Problem 2 is a cutting-plane method, and it worked fine and fast. Reconcile these.

# Add your solution here

5c. Across this whole assignment, exactly one line of defensive code would have caught the FP error in Problem 2 and none of the errors in Problem 4. Name it, and say what a working engineer should conclude about how far solver status checking goes.

# Add your solution here

Declarations¶

Collaboration. List everyone you discussed this assignment with.

AI use. Per the course AI policy, state which tools you used, for which problems, and what you did to verify their output.