Deterministic Global Optimization

Every algorithm in this course so far returns a local solution. When the problem is convex that is enough, because then every local solution is global. When it is not, a local solver answers a question you did not ask — and says nothing about how wrong it might be.

This notebook is the computational companion to the Deterministic Global Optimization handout. The handout derives the method: a convex underestimator gives a lower bound $l$ over a region, a local solve gives an upper bound $u$ , and the region is branched or fathomed according to the gap $u - l$ . Here we compute those bounds, plot them, and watch a global solver do on a real process model what the handout does by hand on a one-dimensional example.

Three parts:

The handout’s example, $\min\ x/4 + \sin x$ on $[-3, 6]$ — the underestimator and the three iterations of spatial branch and bound, computed rather than read off a page.
McCormick envelopes — the relaxation of a bilinear term $w = xy$ , plotted, with the gap formula from the handout verified numerically.
Milk pooling — a bilinear process model with three local optima, where the envelopes of part 2 give a rigorous bound and a global solver closes it.

# 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 numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyomo.environ as pyo

# House figure style. On Colab only the notebook is present, so fall back to the raw URL.
STYLE = (
    "https://raw.githubusercontent.com/ndcbe/optimization/main/figures/dowling.mplstyle"
    if "google.colab" in sys.modules
    else "../../figures/dowling.mplstyle"
)
plt.style.use(STYLE)

# Seeded once, here, for the multistart study in part 3.
rng = np.random.default_rng(seed=0)

solver_lp = pyo.SolverFactory("appsi_highs")
solver_nlp = pyo.SolverFactory("ipopt")

1. The handout’s example¶

\min_x \ f(x) = \frac{x}{4} + \sin x, \qquad -3 \le x \le 6, \qquad \varepsilon = 0.2 .

(1)

The handout builds a convex underestimator out of three straight lines: the global bound $\sin x \ge -1$ , and the tangent lines to $f$ at the two endpoints $x = -3$ and $x = 6$ . Their upper envelope

f^{c}(x) = \max\left\{ \tfrac{1}{4}x - 1, \ \ -0.74x - 3.11, \ \ 1.21x - 6.04 \right\}

(2)

is a maximum of affine functions and therefore convex.

A tangent line to a nonconvex function is not automatically an underestimator, so the handout is careful to say this construction happens to work here. That claim is checkable, and the first thing we do is check it.

def f(x):
    """Objective of the handout's example."""
    return x / 4 + np.sin(x)


def f_prime(x):
    """First derivative of f."""
    return 0.25 + np.cos(x)


def tangent(x0):
    """Slope and intercept of the tangent line to f at x0."""
    slope = f_prime(x0)
    return slope, f(x0) - slope * x0


# The three affine pieces: the sin(x) >= -1 bound, and the two endpoint tangents.
PIECES = [(0.25, -1.0), tangent(-3.0), tangent(6.0)]

for slope, intercept in PIECES:
    print(f"f(x) >= {slope:8.5f} x + {intercept:8.5f}")


def f_under(x, pieces=PIECES):
    """Convex underestimator: the upper envelope of the affine pieces."""
    x = np.asarray(x, dtype=float)
    return np.max([slope * x + intercept for slope, intercept in pieces], axis=0)


# Is it really an underestimator on [-3, 6]? Check on a fine grid.
grid = np.linspace(-3, 6, 200_001)
worst = np.max(f_under(grid) - f(grid))
print(
    f"\nmax(f_under - f) over [-3, 6] = {worst:.3e}  (must be <= 0 to within rounding)"
)
assert worst <= 1e-9, "the piecewise-linear function is NOT an underestimator"

f(x) >=  0.25000 x + -1.00000
f(x) >= -0.73999 x + -3.11110
f(x) >=  1.21017 x + -6.04044

max(f_under - f) over [-3, 6] = 0.000e+00  (must be <= 0 to within rounding)

The maximum of $f^{c} - f$ is zero, attained only at the points of tangency: the relaxation touches $f$ at $x = -3$ and $x = 6$ and lies strictly below it in between. That is the property every bound in this notebook rests on.

The handout’s printed coefficients are rounded to two decimals; the exact tangents are $-0.73999x - 3.11110$ and $1.21017x - 6.04044$ . We use the exact ones so the two lower pieces meet cleanly.

fig, ax = plt.subplots(figsize=(6.4, 4))

x_plot = np.linspace(-3, 6, 400)
ax.plot(x_plot, f(x_plot), label=r"$f(x) = x/4 + \sin x$")
ax.plot(x_plot, f_under(x_plot), label=r"underestimator $f^c(x)$")

ax.axvline(-3, color="0.6", lw=1)
ax.axvline(6, color="0.6", lw=1)
ax.set_xlabel("$x$")
ax.set_ylabel("$f(x)$")
ax.set_xlim(-3.4, 6.4)
ax.legend(loc="upper left")
plt.show()

The three iterations¶

Two problems are solved on every region:

the relaxation $\min f^{c}(x)$ , which is convex — and because $f^c$ is a maximum of affine functions, it is an LP in the epigraph variable $t$ , solved to global optimality by HiGHS;
the original problem $\min f(x)$ , solved locally by Ipopt from a starting point inside the region.

The first gives the lower bound $l$ , the second the upper bound $u$ .

def solve_relaxation(lo, hi):
    """Global minimum of the convex underestimator on [lo, hi], as an LP."""
    m = pyo.ConcreteModel("underestimator LP")
    m.x = pyo.Var(bounds=(lo, hi))
    m.t = pyo.Var(domain=pyo.Reals)
    m.PIECES = pyo.RangeSet(0, len(PIECES) - 1)

    @m.Constraint(m.PIECES)
    def epigraph(m, i):
        slope, intercept = PIECES[i]
        return m.t >= slope * m.x + intercept

    @m.Objective(sense=pyo.minimize)
    def lower_bound(m):
        return m.t

    results = solver_lp.solve(m, load_solutions=False)
    assert pyo.check_optimal_termination(
        results
    ), f"relaxation failed: {results.solver.termination_condition}"
    m.solutions.load_from(results)
    return pyo.value(m.x), pyo.value(m.t)


def solve_local(lo, hi, x0):
    """Local minimum of f on [lo, hi], started from x0."""
    m = pyo.ConcreteModel("original problem")
    m.x = pyo.Var(bounds=(lo, hi), initialize=x0)

    @m.Objective(sense=pyo.minimize)
    def obj(m):
        return m.x / 4 + pyo.sin(m.x)

    results = solver_nlp.solve(m)
    assert pyo.check_optimal_termination(
        results
    ), f"local solve failed: {results.solver.termination_condition}"
    return pyo.value(m.x), pyo.value(m.obj)

EPS = 0.2  # gap tolerance
BRANCH = 4.45971  # the branch point found in iteration 1

regions = [
    ("1: [-3, 6]", -3.0, 6.0, 5.0),
    ("2: A = [-3, 4.46]", -3.0, BRANCH, 0.7),
    ("3: B = [4.46, 6]", BRANCH, 6.0, 5.2),
]

rows = []
for name, lo, hi, x0 in regions:
    x_lower, l = solve_relaxation(lo, hi)
    x_upper, u = solve_local(lo, hi, x0)
    rows.append(
        {
            "iteration / region": name,
            "x^L": x_lower,
            "l": l,
            "x^U": x_upper,
            "u": u,
            "u - l": u - l,
            "decision": "branch" if u - l > EPS else "accept",
        }
    )

iterations = pd.DataFrame(rows).set_index("iteration / region")
display(iterations.round(4))

Loading...

Read the table against the handout.

Iteration 1 finds $\bar{x} = -2.13$ , $l = -1.53$ from the relaxation and $x^U = 4.46$ , $u = 0.147$ from a local solve started on the right. The gap is far larger than $\varepsilon = 0.2$ , so we branch at the upper-bounding point $x^U = 4.46$ — where the relaxation is worst.
Iteration 2 on subregion $A$ keeps the same lower bound (the underestimator did not change) and finds the local solution $x^U = -1.823$ with $u = -1.424$ . The gap is $0.109 < \varepsilon$ : accept.
Iteration 3 on subregion $B$ gives $l = 0.115$ and $u = 0.147$ , a gap of 0.032: accept.

Then $l_B = 0.115 > U = -1.424$ , so subregion $B$ is fathomed — discarded on the strength of its lower bound alone, without ever locating its optimum precisely. The answer is $x^* = -1.823$ in subregion $A$ , and it is global to within $\varepsilon$ .

fig, ax = plt.subplots(figsize=(6.4, 4))

ax.plot(x_plot, f(x_plot), label=r"$f(x)$")
ax.plot(x_plot, f_under(x_plot), label=r"$f^c(x)$")

# the two subregions and the accepted solutions
ax.axvline(BRANCH, color="0.4", lw=1.2)
for name, row in iterations.iloc[1:].iterrows():
    ax.plot(row["x^U"], row["u"], marker="o", color="black", ms=9, ls="none")
    ax.annotate(
        f"$u = {row['u']:.3f}$",
        xy=(row["x^U"], row["u"]),
        xytext=(row["x^U"] + 0.3, row["u"] + 0.45),
        fontsize=12,
    )

ax.annotate("A", xy=(0.5, 1.55), fontsize=16, ha="center")
ax.annotate("B", xy=(5.2, 1.55), fontsize=16, ha="center")
ax.set_xlabel("$x$")
ax.set_ylabel("$f(x)$")
ax.set_xlim(-3.4, 6.4)
ax.set_ylim(-2.2, 2.0)
ax.legend(loc="lower right")
plt.show()

2. McCormick envelopes¶

The one-dimensional example got its relaxation almost for free, because $\sin x \ge -1$ is a bound anyone can write down. Real process models are nonconvex for structural reasons, and the commonest of those is a bilinear term: a component flow is a total flow times a mole fraction, a duty is a flow times a temperature difference.

The handout derives the McCormick envelopes of $w = xy$ on the box $x \in [x^L, x^U]$ , $y \in [y^L, y^U]$ from four sign-definite products:

\begin{aligned} w &\ge x^{L} y + x\, y^{L} - x^{L} y^{L}, &\qquad w &\ge x^{U} y + x\, y^{U} - x^{U} y^{U}, \\ w &\le x^{U} y + x\, y^{L} - x^{U} y^{L}, &\qquad w &\le x^{L} y + x\, y^{U} - x^{L} y^{U}. \end{aligned}

(3)

The first two are the lower envelope, the second two the upper envelope. Below we build them, check that they really do bracket $xy$ , and plot how far apart they are.

def mccormick_envelopes(X, Y, xl, xu, yl, yu):
    """Lower and upper McCormick envelopes of w = x*y on the box, evaluated on a mesh."""
    lower = np.maximum(xl * Y + X * yl - xl * yl, xu * Y + X * yu - xu * yu)
    upper = np.minimum(xu * Y + X * yl - xu * yl, xl * Y + X * yu - xl * yu)
    return lower, upper


xl, xu, yl, yu = 0.0, 1.0, 0.0, 1.0
gx = np.linspace(xl, xu, 41)
gy = np.linspace(yl, yu, 41)
X, Y = np.meshgrid(gx, gy)
W = X * Y
LO, HI = mccormick_envelopes(X, Y, xl, xu, yl, yu)

print(f"lower envelope violated by at most {np.max(LO - W):.2e}")
print(f"upper envelope violated by at most {np.max(W - HI):.2e}")
print(f"maximum separation HI - LO   = {np.max(HI - LO):.4f}")
print(f"(x^U - x^L)(y^U - y^L) / 2   = {0.5 * (xu - xl) * (yu - yl):.4f}")

lower envelope violated by at most 1.11e-16
upper envelope violated by at most 0.00e+00
maximum separation HI - LO   = 0.5000
(x^U - x^L)(y^U - y^L) / 2   = 0.5000

# How far apart the two envelopes are, across the box.
# 7 contour levels, not 20: adjacent bands of any continuous colormap collapse to
# the same grey in print, and 7 keeps them measurably apart.
fig, ax = plt.subplots(figsize=(4.8, 4))
gap = ax.contourf(X, Y, HI - LO, levels=7, cmap="viridis")
ax.contour(X, Y, HI - LO, levels=7, colors="white", linewidths=0.6)
ax.plot(0.5, 0.5, marker="o", color="white", mec="white", ms=9, ls="none")
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_title("envelope gap, upper $-$ lower", fontsize=14)
ax.set_aspect(1)
fig.colorbar(gap, ax=ax, fraction=0.046, pad=0.04)
plt.show()

# The same picture along the diagonal x = y = t, where the gap is widest.
fig, ax = plt.subplots(figsize=(5.2, 4))
t = np.linspace(0, 1, 201)
lo_t, hi_t = mccormick_envelopes(t, t, xl, xu, yl, yu)

# both envelopes share one colour and are told apart by linestyle and a direct
# label: they are two halves of one object, not two independent series
ax.plot(t, t * t, color="black", ls="-")
ax.plot(t, hi_t, color="#0072B2", ls="-.")
ax.plot(t, lo_t, color="#0072B2", ls="--")

ax.annotate("$w = xy$", xy=(0.80, 0.50), fontsize=13)
ax.annotate("upper envelope", xy=(0.06, 0.42), fontsize=12)
ax.annotate("lower envelope", xy=(0.52, 0.05), fontsize=12)
ax.set_xlabel("$x = y = t$")
ax.set_ylabel("$w$")
plt.show()

Two properties matter, and both are visible above. The envelopes are exact on the boundary of the box — the corners the handout mentions, and the edges as well — and loosest at the centre: at $x = y = \tfrac{1}{2}$ they give $0 \le w \le \tfrac{1}{2}$ around a true value of $\tfrac{1}{4}$ . The second figure is the same statement along the diagonal, where the true bilinear function is squeezed between two straight lines that touch it only at the ends.

The handout states that the widest separation between the envelopes is $\tfrac{1}{2}(x^{U} - x^{L})(y^{U} - y^{L})$ — proportional to the product of the box widths. That is why spatial branching works, and it is worth confirming rather than believing.

boxes = [(0, 1, 0, 1), (0, 2, 0, 1), (-1, 3, 2, 5), (0, 0.5, 0, 0.5)]

rows = []
for xl_, xu_, yl_, yu_ in boxes:
    Xb, Yb = np.meshgrid(np.linspace(xl_, xu_, 401), np.linspace(yl_, yu_, 401))
    lo_b, hi_b = mccormick_envelopes(Xb, Yb, xl_, xu_, yl_, yu_)
    rows.append(
        {
            "box": f"[{xl_}, {xu_}] x [{yl_}, {yu_}]",
            "measured max gap": np.max(hi_b - lo_b),
            "dx * dy / 2": 0.5 * (xu_ - xl_) * (yu_ - yl_),
        }
    )

display(pd.DataFrame(rows).set_index("box").round(6))

Loading...

Exact in every case. Halving one side of the box halves the gap; halve both and it quarters. Each spatial split makes every relaxation below it strictly tighter, the lower bounds rise, and regions start failing the fathoming test — which is the entire economy of spatial branch and bound.

3. Milk pooling¶

Now a process model. A bulk distributor blends raw milk to meet customer specifications on fat content. Milk from two local farms can be delivered separately, but the two remote farms share a single truck with one tank, so their milk must be mixed into a pool of uniform composition $p$ before it reaches the blending station.

That single sentence is what makes the problem nonconvex. The fat delivered to customer $c$ from the pool is $p\,y_c$ — pool composition times pool flow — a bilinear term, one intensive variable times one extensive variable. It is the same structure as $w = xy$ in part 2, and it is the reason pooling problems have been a benchmark for global optimization for nearly fifty years.

customers = pd.DataFrame(
    {
        "Customer 1": {"min_fat": 0.045, "price": 52.0, "demand": 6000.0},
        "Customer 2": {"min_fat": 0.030, "price": 48.0, "demand": 2500.0},
        "Customer 3": {"min_fat": 0.040, "price": 50.0, "demand": 4000.0},
    }
).T

suppliers = pd.DataFrame(
    {
        "Farm A": {"fat": 0.045, "cost": 45.0, "location": "local"},
        "Farm B": {"fat": 0.030, "cost": 42.0, "location": "local"},
        "Farm C": {"fat": 0.033, "cost": 37.0, "location": "remote"},
        "Farm D": {"fat": 0.050, "cost": 45.0, "location": "remote"},
    }
).T

LOCAL = suppliers.index[suppliers["location"] == "local"].tolist()
REMOTE = suppliers.index[suppliers["location"] == "remote"].tolist()

print("Customers")
display(customers)
print("Suppliers")
display(suppliers)

Customers

Loading...

Suppliers

Loading...

Two linear reference cases¶

Before the nonconvex model, two LPs that bracket it.

Local farms only — business as usual. No pool, no bilinear term.
Every farm, no pooling — what the distributor could earn with a second truck, so that remote milk arrives unmixed. This relaxes the pooling restriction, so it is an upper bound on the pooling profit.

Both are ordinary blending LPs: buy from each source, deliver to each customer, respect demand and a minimum fat fraction.

def build_blending_lp(sources):
    """Blending LP over a given set of supplier names.

    Arguments:
        sources: list of supplier names that can deliver to customers separately

    Returns:
        a Pyomo ConcreteModel
    """
    m = pyo.ConcreteModel("milk blending LP")
    m.S = pyo.Set(initialize=sources)
    m.C = pyo.Set(initialize=customers.index.tolist())

    # z[s, c] = milk shipped from supplier s into the blend for customer c
    m.z = pyo.Var(m.S, m.C, domain=pyo.NonNegativeReals)

    @m.Constraint(m.C)
    def demand(m, c):
        return sum(m.z[s, c] for s in m.S) <= customers.loc[c, "demand"]

    @m.Constraint(m.C)
    def fat_spec(m, c):
        # linear blending: delivered fat >= required fraction * delivered volume
        return sum(suppliers.loc[s, "fat"] * m.z[s, c] for s in m.S) >= customers.loc[
            c, "min_fat"
        ] * sum(m.z[s, c] for s in m.S)

    @m.Objective(sense=pyo.maximize)
    def profit(m):
        return sum(
            m.z[s, c] * (customers.loc[c, "price"] - suppliers.loc[s, "cost"])
            for s in m.S
            for c in m.C
        )

    return m


reference = {}
for label, sources in [
    ("local farms only", LOCAL),
    ("all farms, no pooling", suppliers.index.tolist()),
]:
    m = build_blending_lp(sources)
    results = solver_lp.solve(m, load_solutions=False)
    assert pyo.check_optimal_termination(
        results
    ), f"{label}: {results.solver.termination_condition}"
    m.solutions.load_from(results)
    reference[label] = pyo.value(m.profit)
    print(f"{label:>24s}: profit = {reference[label]:10,.0f}")

        local farms only: profit =     81,000

   all farms, no pooling: profit =    122,441

The pooling model¶

With the pool, let $x_r$ be the milk bought from remote farm $r$ , $y_c$ the pooled milk delivered to customer $c$ , $z_{l,c}$ the local milk delivered to customer $c$ , and $p$ the fat fraction of the pool:

\begin{aligned} \max \quad & \sum_{l,c} (\text{price}_c - \text{cost}_l) z_{l,c} + \sum_c \text{price}_c\, y_c - \sum_r \text{cost}_r\, x_r \\ \text{s.t.}\quad & \sum_l z_{l,c} + y_c \le \text{demand}_c && \forall c \\ & \sum_r x_r = \sum_c y_c \\ & \sum_r \text{fat}_r\, x_r = \underbrace{p \sum_c y_c}_{\text{bilinear}} \\ & \underbrace{p\, y_c}_{\text{bilinear}} + \sum_l \text{fat}_l\, z_{l,c} \ \ge\ \text{fat}^{\min}_c \Big( \sum_l z_{l,c} + y_c \Big) && \forall c \\ & p, x_r, y_c, z_{l,c} \ge 0 . \end{aligned}

(4)

Two bilinear terms, both products of the pool composition with a flow. Everything else is linear.

The build function below takes p either as a fixed Param — which makes the model an LP — or as a Var, which makes it a nonconvex NLP. Same constraints either way, which is what lets us scan $p$ and then optimize over it.

P_LO = suppliers.loc[REMOTE, "fat"].min()
P_HI = suppliers.loc[REMOTE, "fat"].max()


def build_pooling(p=None, p_init=None):
    """Milk pooling model in the p-parameterization.

    Arguments:
        p: pool fat fraction. If a float, it is fixed as a Param and the model is an LP.
           If None, p becomes a decision variable and the model is a bilinear NLP.
        p_init: initial value for the p variable (ignored when p is fixed)

    Returns:
        a Pyomo ConcreteModel
    """
    m = pyo.ConcreteModel("milk pooling")
    m.L = pyo.Set(initialize=LOCAL)
    m.R = pyo.Set(initialize=REMOTE)
    m.C = pyo.Set(initialize=customers.index.tolist())

    m.x = pyo.Var(m.R, domain=pyo.NonNegativeReals)  # bought from remote farm r
    m.y = pyo.Var(m.C, bounds=lambda m, c: (0, customers.loc[c, "demand"]))
    m.z = pyo.Var(m.L, m.C, domain=pyo.NonNegativeReals)

    if p is None:
        m.p = pyo.Var(
            bounds=(P_LO, P_HI),
            initialize=0.5 * (P_LO + P_HI) if p_init is None else p_init,
        )
    else:
        m.p = pyo.Param(initialize=p)

    @m.Constraint(m.C)
    def demand(m, c):
        return sum(m.z[l, c] for l in m.L) + m.y[c] <= customers.loc[c, "demand"]

    @m.Constraint()
    def pool_balance(m):
        return sum(m.x[r] for r in m.R) == sum(m.y[c] for c in m.C)

    @m.Constraint()
    def pool_fat(m):
        # fat into the pool = fat out of the pool             [bilinear]
        return sum(suppliers.loc[r, "fat"] * m.x[r] for r in m.R) == m.p * sum(
            m.y[c] for c in m.C
        )

    @m.Constraint(m.C)
    def fat_spec(m, c):
        # fat delivered to customer c meets the specification [bilinear]
        return m.p * m.y[c] + sum(
            suppliers.loc[l, "fat"] * m.z[l, c] for l in m.L
        ) >= customers.loc[c, "min_fat"] * (sum(m.z[l, c] for l in m.L) + m.y[c])

    @m.Objective(sense=pyo.maximize)
    def profit(m):
        return (
            sum(
                m.z[l, c] * (customers.loc[c, "price"] - suppliers.loc[l, "cost"])
                for l in m.L
                for c in m.C
            )
            + sum(customers.loc[c, "price"] * m.y[c] for c in m.C)
            - sum(suppliers.loc[r, "cost"] * m.x[r] for r in m.R)
        )

    return m

What the nonconvexity looks like¶

Fix $p$ and the model is an LP. So sweep $p$ across its range, solve an LP at each value, and plot the profit. This is not an algorithm — it only works because there is exactly one troublesome variable — but it shows what the solver is up against.

p_scan = np.linspace(P_LO, P_HI, 201)
profit_scan = np.full_like(p_scan, np.nan)

for i, p in enumerate(p_scan):
    m = build_pooling(p=p)
    results = solver_lp.solve(m, load_solutions=False)
    if pyo.check_optimal_termination(results):
        m.solutions.load_from(results)
        profit_scan[i] = pyo.value(m.profit)
    else:
        print(f"p = {p:.4f}: {results.solver.termination_condition}")


def local_maxima(values):
    """Indices of local maxima of a 1-D array, endpoints included."""
    out = []
    for i in range(len(values)):
        left = values[i - 1] if i > 0 else -np.inf
        right = values[i + 1] if i < len(values) - 1 else -np.inf
        if values[i] >= left and values[i] >= right and values[i] > min(left, right):
            out.append(i)
    return out


maxima = local_maxima(profit_scan)
print("local maxima found by the scan:")
for i in maxima:
    print(f"  p = {p_scan[i]:.4f}   profit = {profit_scan[i]:10,.0f}")

local maxima found by the scan:
  p = 0.0330   profit =    102,833
  p = 0.0400   profit =    100,067
  p = 0.0451   profit =    101,134

fig, ax = plt.subplots(figsize=(6.4, 4))

ax.plot(p_scan, profit_scan / 1000, label="profit at fixed $p$")
ax.plot(
    p_scan[maxima],
    profit_scan[maxima] / 1000,
    marker="o",
    ls="none",
    color="black",
    ms=9,
    label="local maxima",
)

ax.set_xlabel("pool fat fraction $p$")
ax.set_ylabel("profit [thousand \\$]")
ax.legend(loc="lower left", fontsize=11)
plt.show()

Three local maxima, separated by two local minima. A local NLP solver will converge to one of them and report success. Which one depends entirely on where it starts.

A rigorous bound from McCormick envelopes¶

Replace each bilinear term $p\,y_c$ by a new variable $w_c$ and add the four envelope inequalities from part 2, with $x \to p$ and $y \to y_c$ . Every remaining constraint is linear, so the relaxation is an LP — solved to global optimality — and because it relaxes the original problem, its optimal value is a valid upper bound on the pooling profit.

The bounds the envelopes need are already in the model: $p \in [0.033, 0.050]$ from the remote farms’ fat contents, and $y_c \in [0, \text{demand}_c]$ .

def build_pooling_mccormick():
    """McCormick relaxation of the pooling model: an LP upper bound on the profit."""
    m = build_pooling(p=None)

    # w[c] replaces the bilinear product p * y[c]
    m.w = pyo.Var(m.C, domain=pyo.NonNegativeReals)

    # the two constraints containing bilinear terms are rebuilt in terms of w
    m.del_component(m.pool_fat)
    m.del_component(m.fat_spec)

    @m.Constraint()
    def pool_fat(m):
        return sum(suppliers.loc[r, "fat"] * m.x[r] for r in m.R) == sum(
            m.w[c] for c in m.C
        )

    @m.Constraint(m.C)
    def fat_spec(m, c):
        return m.w[c] + sum(
            suppliers.loc[l, "fat"] * m.z[l, c] for l in m.L
        ) >= customers.loc[c, "min_fat"] * (sum(m.z[l, c] for l in m.L) + m.y[c])

    # McCormick envelopes of w[c] = p * y[c] on [p_lo, p_hi] x [0, demand_c]
    @m.Constraint(m.C)
    def mccormick_lower_1(m, c):
        return m.w[c] >= m.p.lb * m.y[c] + m.y[c].lb * m.p - m.p.lb * m.y[c].lb

    @m.Constraint(m.C)
    def mccormick_lower_2(m, c):
        return m.w[c] >= m.p.ub * m.y[c] + m.y[c].ub * m.p - m.p.ub * m.y[c].ub

    @m.Constraint(m.C)
    def mccormick_upper_1(m, c):
        return m.w[c] <= m.p.ub * m.y[c] + m.y[c].lb * m.p - m.p.ub * m.y[c].lb

    @m.Constraint(m.C)
    def mccormick_upper_2(m, c):
        return m.w[c] <= m.p.lb * m.y[c] + m.y[c].ub * m.p - m.p.lb * m.y[c].ub

    return m


m_relaxed = build_pooling_mccormick()
results = solver_lp.solve(m_relaxed, load_solutions=False)
assert pyo.check_optimal_termination(
    results
), f"relaxation failed: {results.solver.termination_condition}"
m_relaxed.solutions.load_from(results)

relaxed_profit = pyo.value(m_relaxed.profit)
relaxed_p = pyo.value(m_relaxed.p)
print(f"McCormick upper bound = {relaxed_profit:10,.0f}  at p = {relaxed_p:.4f}")
print(f"best profit seen in the scan = {np.nanmax(profit_scan):10,.0f}")
print(f"remaining gap = {relaxed_profit - np.nanmax(profit_scan):,.0f}")

McCormick upper bound =    111,412  at p = 0.0400
best profit seen in the scan =    102,833
remaining gap = 8,578

The relaxation returns 111,412 — a number no feasible operating point can beat. It is a genuine certificate, and it is loose: the best profit actually achievable is around 103,000. The relaxation also reports $p = 0.040$ , which is not the composition of any of the three local maxima.

That combination is exactly what the handout warns about, and it is worth being blunt: a relaxation is a bound, not a plan. Taking $p = 0.040$ from the relaxation and running the true model there is a legitimate heuristic, but it lands on the worst of the three local maxima.

Closing that gap is what spatial branch and bound does: split the range of $p$ , rebuild the envelopes on each half — the gap shrinks in proportion to the box width, as measured in part 2 — and repeat until the bound meets an achievable profit.

# What does the relaxation's p actually deliver?
m_at_relaxed_p = build_pooling(p=relaxed_p)
results = solver_lp.solve(m_at_relaxed_p, load_solutions=False)
assert pyo.check_optimal_termination(results), results.solver.termination_condition
m_at_relaxed_p.solutions.load_from(results)
achievable_at_relaxed_p = pyo.value(m_at_relaxed_p.profit)
print(
    f"profit actually achievable at p = {relaxed_p:.4f}: {achievable_at_relaxed_p:,.0f}"
)

profit actually achievable at p = 0.0400: 100,088

Local solvers, and how much the starting point matters¶

Ipopt solves the bilinear model to a point satisfying the KKT conditions. Run it from twelve random starting values of $p$ — the random number generator was seeded at the top of the notebook, so this is reproducible — and count where it lands.

multistart = []
for trial in range(12):
    p0 = rng.uniform(P_LO, P_HI)
    m = build_pooling(p=None, p_init=p0)
    results = solver_nlp.solve(m)
    tc = results.solver.termination_condition
    if pyo.check_optimal_termination(results):
        multistart.append(
            {"p start": p0, "p found": pyo.value(m.p), "profit": pyo.value(m.profit)}
        )
    else:
        multistart.append({"p start": p0, "p found": np.nan, "profit": np.nan})
        print(f"trial {trial}: {tc}")

multistart = pd.DataFrame(multistart)
display(multistart.round(4))

print("\ndistinct local optima reached:")
display(
    multistart.round(4)
    .groupby("profit")
    .agg(**{"p found": ("p found", "first"), "count": ("p found", "size")})
)

Loading...


distinct local optima reached:

Loading...

Three distinct answers, from twelve runs of the same solver on the same model. Every one of them terminated optimally; every one of them is a local maximum; and Ipopt has no way to tell you which is which. The termination condition tells you the solve succeeded. It does not tell you the answer is global — and no local optimality condition ever can, because optimality conditions are statements about derivatives at a point.

A global solver¶

Couenne is a deterministic global solver: it does what the handout’s flowchart does — relax, bound, branch, fathom — on a spatial branch and bound tree, with the McCormick envelopes of part 2 among its relaxations. It returns the global optimum and a proof that nothing better exists.

m_global = build_pooling(p=None)
results = pyo.SolverFactory("couenne").solve(m_global)
assert pyo.check_optimal_termination(
    results
), f"couenne: {results.solver.termination_condition}"

global_profit = pyo.value(m_global.profit)
print(
    f"global optimum: profit = {global_profit:,.2f} at p = {pyo.value(m_global.p):.4f}"
)

summary = pd.DataFrame(
    [
        ("local farms only (LP)", reference["local farms only"], "achievable"),
        ("global optimum with pooling (Couenne)", global_profit, "achievable"),
        ("best local optimum found by Ipopt", multistart["profit"].max(), "achievable"),
        (
            "worst local optimum found by Ipopt",
            multistart["profit"].min(),
            "achievable",
        ),
        ("McCormick relaxation (LP)", relaxed_profit, "upper bound"),
        (
            "all farms, no pooling (LP)",
            reference["all farms, no pooling"],
            "upper bound",
        ),
    ],
    columns=["case", "profit", "kind"],
).set_index("case")

display(summary.round(1))

global optimum: profit = 102,833.33 at p = 0.0330

Loading...

fig, ax = plt.subplots(figsize=(6.4, 4))

ax.plot(p_scan, profit_scan / 1000, label="profit at fixed $p$")
ax.axhline(relaxed_profit / 1000, color="0.35", ls="--", lw=2, label="McCormick bound")
ax.plot(
    pyo.value(m_global.p),
    global_profit / 1000,
    marker="o",
    ls="none",
    color="black",
    ms=11,
    label="global optimum",
)
ax.plot(
    multistart["p found"],
    multistart["profit"] / 1000,
    marker="s",
    ls="none",
    mfc="none",
    mec="black",
    ms=10,
    label="Ipopt, 12 random starts",
)

ax.set_xlabel("pool fat fraction $p$")
ax.set_ylabel("profit [thousand \\$]")
ax.set_ylim(85, 114)
ax.legend(loc="lower left", fontsize=10)
plt.show()

The picture is the whole lecture in one axes. The solid curve is the true profit as a function of $p$ , with three local maxima. The dashed line is a rigorous upper bound obtained from an LP, valid over the entire range without evaluating the true model anywhere. The open squares are where a local solver lands from twelve different starting points. The filled circle is the global optimum, and Couenne can prove it.

Three closing points, all of which the handout makes and this notebook only illustrates:

Pooling costs money. With a second truck the distributor would earn 122,441; pooling caps the profit at 102,833. That difference — about 19,600 a year — is what the second truck is worth, and the comparison is only meaningful because 102,833 is a proved optimum rather than whatever a local solver happened to find.
A relaxation bounds, it does not plan. The McCormick bound of 111,412 is valid and unachievable, and the composition it reports leads to the worst of the three local maxima.
Global costs time. Couenne solves this in seconds; on a full refinery pooling model it would not. The tolerance $\varepsilon$ is where you buy back that time, and an answer proved to be within $\varepsilon$ of global is a categorically different object from an answer that merely converged.

1. The handout’s example¶

The three iterations¶

2. McCormick envelopes¶

3. Milk pooling¶

Two linear reference cases¶

The pooling model¶

What the nonconvexity looks like¶

A rigorous bound from McCormick envelopes¶

Local solvers, and how much the starting point matters¶

A global solver¶

Further reading¶