Milk Pooling and Blending

Adapted, with thanks, from work by the late Jeffrey Kantor.

This notebook is adapted from notebook 5.1, Milk pooling and blending, in the companion notebook collection for M. Postek, A. Zocca, J. Gromicho and J. Kantor, Hands-On Mathematical Optimization with Python (Cambridge University Press, 2025). Jeff Kantor was a colleague in this department; the example, the data and the framing of the three options are his.

The underlying problem is the canonical pooling benchmark of C. A. Haverly, “Studies of the behavior of recursion for the pooling problem,” ACM SIGMAP Bulletin 25 (1978), 19--28.

MIT License

Copyright (c) 2022 Jeffrey Kantor

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

What is different here. The models are rewritten with named pyo.Set and pyo.Param components so that the Pyomo code reads as the set-notation model in the lecture notes, the solver is ipopt throughout rather than appsi_highs, and the convex (McCormick) relaxation is left to the global optimization lecture. Two typographical errors in the source notebook’s mathematics are corrected: its blending constraints are printed with the summation running over the cross product $L \times C$ inside a constraint that is already quantified for all $c \in C$ , and the sum must run over $l \in L$ alone. Its Pyomo code is correct as written.

# 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
from pyomo.environ import units as u
from pyomo.util.check_units import assert_units_consistent

# The source data is scale-free: milk is measured in "units" and priced per unit.
# Pyomo's unit library has neither, so declare both. Fat content stays
# dimensionless -- it is a mass fraction, kg of fat per kg of milk.
u.load_definitions_from_strings(["USD = [currency]", "milk = [milk]"])

Problem description¶

A bulk distributor supplies custom milk blends to several customers. Each customer specifies a minimum fat content, pays a fixed price per unit, and will buy no more than a stated maximum. The distributor sources raw milk from local farms, each producing milk of known fat content at a known cost.

The distributor has found cheaper raw milk at some remote farms. But it owns only one truck with a single tank for the remote route, so milk from the remote farms must be mixed in that tank before transport. That creates a pool of uniform composition, which is then blended with local milk to meet each customer’s requirement.

Three options:

Business as usual --- local farms only.
Buy a second truck --- keep each remote farm’s milk separate, so everything blends freely.
Pool the remote farms --- one truck, one tank, one composition.

Options 1 and 2 are linear programs. Option 3 is not, and that is the point of the example.

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

# Farm and customer names
local_farms = suppliers.index[suppliers["location"] == "local"].tolist()
remote_farms = suppliers.index[suppliers["location"] == "remote"].tolist()
all_farms = local_farms + remote_farms
customer_names = customers.index.tolist()

# Problem data as plain dictionaries, keyed by farm or customer name

# Fat fraction of the milk from farm s, phi_s [dimensionless]
fat = suppliers["fat"].astype(float).to_dict()

# Purchase cost per unit from farm s, kappa_s [USD/milk]
cost = suppliers["cost"].astype(float).to_dict()

# Price paid per unit by customer k, pi_k [USD/milk]
price = customers["price"].astype(float).to_dict()

# Maximum demand of customer k, d_k [milk]
demand = customers["demand"].astype(float).to_dict()

# Minimum fat fraction customer k accepts, phi_min_k [dimensionless]
min_fat = customers["min_fat"].astype(float).to_dict()

print("Customers")
display(customers)

print("Suppliers")
display(suppliers)

Customers

Suppliers

Options 1 and 2 --- blending only, no pool¶

Let $\mathcal{S}$ be a set of farms and $\mathcal{K}$ the customers. The decision variable $z_{s,k}$ is the amount of milk from farm $s$ blended into the delivery to customer $k$ .

\begin{align} \max_{z} \quad & \sum_{s \in \mathcal{S}} \sum_{k \in \mathcal{K}} (\pi_k - \kappa_s) z_{s,k} && \text{profit} \\ \text{s.t.} \quad & \sum_{s \in \mathcal{S}} z_{s,k} \leq d_k && \forall k \in \mathcal{K} \quad \text{demand} \\ & \sum_{s \in \mathcal{S}} \phi_s z_{s,k} \geq \phi^{\min}_k \sum_{s \in \mathcal{S}} z_{s,k} && \forall k \in \mathcal{K} \quad \text{blend quality} \\ & z_{s,k} \geq 0 \end{align}

(1)

Here $\phi_s$ is the fat fraction of farm $s$ ’s milk, $\kappa_s$ its cost per unit, $\pi_k$ the price customer $k$ pays, $d_k$ its maximum demand and $\phi^{\min}_k$ its minimum acceptable fat fraction. Option 1 takes $\mathcal{S} = \mathcal{L}$ , the local farms. Option 2 takes $\mathcal{S} = \mathcal{L} \cup \mathcal{R}$ , all of them.

Every term is linear, so this is a linear program.

Units. Milk here is measured in generic units, so the notebook declares a milk unit and a USD currency and hangs them on every Param and Var. Fat content stays dimensionless, because it is a mass fraction. Both builders end with assert_units_consistent, which is the part that earns the declaration: it is what would catch a price or a cost mixed into a constraint that balances amounts of milk. It cannot catch everything: fat content is dimensionless, so a fat balance and a milk balance carry the same units.

def create_blending_model(farms):
    """Create the linear blending model (no pool) in Pyomo

    Arguments:
        farms: list of farm names to source from

    Returns:
        m: Pyomo concrete model
    """

    m = pyo.ConcreteModel("Milk blending")

    # Farms, S in the notes
    m.S = pyo.Set(initialize=farms)

    # Customers, K in the notes
    m.K = pyo.Set(initialize=customer_names)

    # Fat fraction of the milk from farm s, phi_s [dimensionless]
    m.fat = pyo.Param(m.S, initialize={s: fat[s] for s in farms})

    # Purchase cost per unit from farm s, kappa_s [USD/milk]
    m.cost = pyo.Param(
        m.S, initialize={s: cost[s] for s in farms}, units=u.USD / u.milk
    )

    # Price paid per unit by customer k, pi_k [USD/milk]
    m.price = pyo.Param(m.K, initialize=price, units=u.USD / u.milk)

    # Maximum demand of customer k, d_k [milk]
    m.demand = pyo.Param(m.K, initialize=demand, units=u.milk)

    # Minimum fat fraction customer k accepts, phi_min_k [dimensionless]
    m.min_fat = pyo.Param(m.K, initialize=min_fat)

    # Milk shipped from farm s to customer k, z_{s,k} [milk]
    m.z = pyo.Var(m.S, m.K, domain=pyo.NonNegativeReals, units=u.milk)

    # Maximize profit [USD]
    @m.Objective(sense=pyo.maximize)
    def profit(b):
        return sum((b.price[k] - b.cost[s]) * b.z[s, k] for s in b.S for k in b.K)

    # Deliveries cannot exceed each customer's demand [milk]
    @m.Constraint(m.K)
    def demand_limit(b, k):
        return sum(b.z[s, k] for s in b.S) <= b.demand[k]

    # Linear blending: each delivery meets its fat spec [milk]
    @m.Constraint(m.K)
    def blend_quality(b, k):
        return sum(b.fat[s] * b.z[s, k] for s in b.S) >= b.min_fat[k] * sum(
            b.z[s, k] for s in b.S
        )

    # Raises UnitsError if any constraint or the objective is inconsistent
    assert_units_consistent(m)

    return m

solver = pyo.SolverFactory("ipopt")


def report_blending(m, title):
    results = solver.solve(m)
    assert pyo.check_optimal_termination(results), f"Solve failed for {title}"
    print(f"{title}: profit = {pyo.value(m.profit):,.2f}\n")

    rows = []
    for k in m.K:
        row = {s: pyo.value(m.z[s, k]) for s in m.S}
        total = sum(row.values())
        row["Total"] = total
        row["fat delivered"] = (
            sum(pyo.value(m.z[s, k]) * pyo.value(m.fat[s]) for s in m.S) / total
            if total > 0
            else np.nan
        )
        row["min_fat"] = pyo.value(m.min_fat[k])
        rows.append(row)
    display(pd.DataFrame(rows, index=list(m.K)).round(4))
    return pyo.value(m.profit)


profit_1 = report_blending(
    create_blending_model(local_farms), "Option 1, local farms only"
)

Option 1, local farms only: profit = 81,000.00

profit_2 = report_blending(
    create_blending_model(all_farms), "Option 2, all farms, one truck each"
)

Option 2, all farms, one truck each: profit = 122,441.18

Sourcing from the remote farms is worth a great deal --- but Option 2 needs a second truck to keep the two remote supplies separated in transit. Note also that the local farms are displaced entirely, and that Customer 2 is delivered milk richer than it requires. That is product giveaway: the optimizer hands over quality because the alternative costs more.

Option 3 --- the pooling problem¶

Only one truck is available for the remote route, so purchases from the remote farms are combined into a single pool of uniform composition, transported, and then blended with local milk.

This is the $p$ -parameterization of the pooling problem: the pool composition is itself a decision variable, $p$ .

Sets. $\mathcal{L}$ local farms; $\mathcal{R}$ remote farms; $\mathcal{K}$ customers.

Variables. $z_{l,k}$ , local milk shipped from $l \in \mathcal{L}$ to customer $k$ ; $x_r$ , milk bought from remote farm $r$ ; $y_k$ , pooled milk delivered to customer $k$ ; and $p$ , the fat fraction of the pool.

\begin{align} \max_{x, y, z, p} \quad & \sum_{l \in \mathcal{L}} \sum_{k \in \mathcal{K}} (\pi_k - \kappa_l) z_{l,k} + \sum_{k \in \mathcal{K}} \pi_k y_k - \sum_{r \in \mathcal{R}} \kappa_r x_r && \text{profit} \\ \text{s.t.} \quad & \sum_{l \in \mathcal{L}} z_{l,k} + y_k \leq d_k && \forall k \in \mathcal{K} \quad \text{demand} \\ & \sum_{r \in \mathcal{R}} x_r = \sum_{k \in \mathcal{K}} y_k && \text{pool balance} \\ & \sum_{r \in \mathcal{R}} \phi_r x_r = \underbrace{p \sum_{r \in \mathcal{R}} x_r}_{\textbf{bilinear}} && \text{pool fat balance} \\ & \underbrace{p\, y_k}_{\textbf{bilinear}} + \sum_{l \in \mathcal{L}} \phi_l z_{l,k} \geq \phi^{\min}_k \Big( \sum_{l \in \mathcal{L}} z_{l,k} + y_k \Big) && \forall k \in \mathcal{K} \quad \text{blend quality} \end{align}

(2)

Every variable is nonnegative, $y_k \leq d_k$ , and $p$ is bounded below and above by the smallest and largest fat fractions among the remote farms --- a pool cannot be richer than its richest input.

The pool fat balance can be written either as $p \sum_{r \in \mathcal{R}} x_r$ , as above and as the Pyomo code below does, or as $p \sum_{k \in \mathcal{K}} y_k$ . The two agree on the feasible set, because the pool balance constraint immediately above forces $\sum_{r} x_r = \sum_{k} y_k$ .

The last two constraints contain bilinear terms: a product of the composition variable $p$ with a flow variable. A bilinear term is smooth and differentiable everywhere, and irreducibly nonconvex. Unlike an absolute value, no change of variables removes it.

def create_pooling_model():
    """Create the bilinear milk pooling model in Pyomo

    Returns:
        m: Pyomo concrete model
    """

    m = pyo.ConcreteModel("Milk pooling")

    # Local farms, L in the notes
    m.L = pyo.Set(initialize=local_farms)

    # Remote farms, R in the notes
    m.R = pyo.Set(initialize=remote_farms)

    # All farms, L union R
    m.S = m.L | m.R

    # Customers, K in the notes
    m.K = pyo.Set(initialize=customer_names)

    # Fat fraction of the milk from farm s, phi_s [dimensionless]
    m.fat = pyo.Param(m.S, initialize=fat)

    # Purchase cost per unit from farm s, kappa_s [USD/milk]
    m.cost = pyo.Param(m.S, initialize=cost, units=u.USD / u.milk)

    # Price paid per unit by customer k, pi_k [USD/milk]
    m.price = pyo.Param(m.K, initialize=price, units=u.USD / u.milk)

    # Maximum demand of customer k, d_k [milk]
    m.demand = pyo.Param(m.K, initialize=demand, units=u.milk)

    # Minimum fat fraction customer k accepts, phi_min_k [dimensionless]
    m.min_fat = pyo.Param(m.K, initialize=min_fat)

    # Local milk shipped from farm l to customer k, z_{l,k} [milk]
    m.z = pyo.Var(m.L, m.K, domain=pyo.NonNegativeReals, units=u.milk)

    # Milk bought from remote farm r, x_r [milk]
    m.x = pyo.Var(m.R, domain=pyo.NonNegativeReals, units=u.milk)

    # Pooled milk delivered to customer k, y_k [milk]
    m.y = pyo.Var(m.K, bounds=lambda m, k: (0, m.demand[k]), units=u.milk)

    # A pool cannot be richer than its richest input, nor leaner than its
    # leanest.
    # Fat fraction of the pool, p [dimensionless]
    m.p = pyo.Var(
        bounds=(min(fat[r] for r in remote_farms), max(fat[r] for r in remote_farms))
    )

    # Maximize profit [USD]
    @m.Objective(sense=pyo.maximize)
    def profit(b):
        return (
            sum((b.price[k] - b.cost[l]) * b.z[l, k] for l in b.L for k in b.K)
            + sum(b.price[k] * b.y[k] for k in b.K)
            - sum(b.cost[r] * b.x[r] for r in b.R)
        )

    # Local milk plus pooled milk cannot exceed each customer's demand [milk]
    @m.Constraint(m.K)
    def demand_limit(b, k):
        return sum(b.z[l, k] for l in b.L) + b.y[k] <= b.demand[k]

    # What is bought from the remote farms is what leaves the pool [milk]
    @m.Constraint()
    def pool_balance(b):
        return sum(b.x[r] for r in b.R) == sum(b.y[k] for k in b.K)

    # Fat into the pool equals fat out. BILINEAR. [milk]
    @m.Constraint()
    def pool_quality(b):
        return sum(b.fat[r] * b.x[r] for r in b.R) == b.p * sum(b.x[r] for r in b.R)

    # Each delivery meets its fat spec. BILINEAR. [milk]
    @m.Constraint(m.K)
    def blend_quality(b, k):
        fat_delivered = b.p * b.y[k] + sum(b.fat[l] * b.z[l, k] for l in b.L)
        milk_delivered = sum(b.z[l, k] for l in b.L) + b.y[k]
        return fat_delivered >= b.min_fat[k] * milk_delivered

    # Raises UnitsError if any constraint or the objective is inconsistent
    assert_units_consistent(m)

    return m

def report_pooling(m, title):
    print(f"{title}")
    print(f"  pool composition p = {pyo.value(m.p):6.4f}")
    print(f"  profit             = {pyo.value(m.profit):,.2f}\n")

    rows = []
    for s in m.S:
        row = {"fat": pyo.value(m.fat[s]), "cost": pyo.value(m.cost[s])}
        for k in m.K:
            row[k] = pyo.value(m.z[s, k]) if s in m.L else 0.0
        row["Pool"] = pyo.value(m.x[s]) if s in m.R else 0.0
        rows.append(row)
    print("  Supplier report")
    display(pd.DataFrame(rows, index=list(m.S)).round(4))

    rows = []
    for k in m.K:
        row = {"min_fat": pyo.value(m.min_fat[k]), "demand": pyo.value(m.demand[k])}
        for l in m.L:
            row[l] = pyo.value(m.z[l, k])
        row["Pool"] = pyo.value(m.y[k])
        total = sum(pyo.value(m.z[l, k]) for l in m.L) + pyo.value(m.y[k])
        row["Total"] = total
        row["fat delivered"] = (
            (
                sum(pyo.value(m.z[l, k]) * pyo.value(m.fat[l]) for l in m.L)
                + pyo.value(m.y[k]) * pyo.value(m.p)
            )
            / total
            if total > 0
            else np.nan
        )
        rows.append(row)
    print("  Customer report")
    display(pd.DataFrame(rows, index=list(m.K)).round(4))

Solve it twice¶

bound_frac is an Ipopt option that controls only how far inside its bounds the solver places the initial point. It changes nothing about the model. Solve the same model with two values of it.

nlp_solver = pyo.SolverFactory("ipopt")

nlp_solver.options["bound_frac"] = 0.01
m_a = create_pooling_model()
results = nlp_solver.solve(m_a)
assert pyo.check_optimal_termination(results), "Solve failed"
report_pooling(m_a, "Option 3, pooled, bound_frac = 0.01")

Option 3, pooled, bound_frac = 0.01
  pool composition p = 0.0330
  profit             = 102,833.33

  Supplier report

  Customer report

nlp_solver.options["bound_frac"] = 0.5
m_b = create_pooling_model()
results = nlp_solver.solve(m_b)
assert pyo.check_optimal_termination(results), "Solve failed"
report_pooling(m_b, "Option 3, pooled, bound_frac = 0.5")

Option 3, pooled, bound_frac = 0.5
  pool composition p = 0.0450
  profit             = 101,392.16

  Supplier report

  Customer report

Activity

Same model, same solver, same tolerances. One option changed, and it only affects where the solver starts. Which of the two answers is right, and how would you know?

Why: fix $p$ and look¶

If $p$ is held fixed, every remaining constraint is linear and the problem is an LP. Sweeping $p$ over its range and solving the LP at each value shows the profit as a function of the pool composition. It has three local maxima.

It is tempting to say that Ipopt converges to whichever one it starts nearest. That is false, and this model is the counterexample: reading Ipopt’s pushed-in initial point back with max_iter = 0, bound_frac = 0.5 starts $p$ at 0.0415 and bound_frac = 0.4 starts it at 0.0398 --- essentially on the local maximum near 0.040 --- and both converge to 0.045. The reason is the teachable part: bound_frac pushes every variable off its bounds, so the twelve-vector the barrier method actually starts from is nowhere near the local solution whose $p$ value it happens to share. Proximity in one coordinate predicts nothing.

def profit_at_fixed_p(p_value):
    """Solve the pooling model with the pool composition p fixed. This is an LP."""
    m = create_pooling_model()
    m.p.fix(p_value)
    results = solver.solve(m)
    if not pyo.check_optimal_termination(results):
        return np.nan
    return pyo.value(m.profit)


p_grid = np.linspace(0.033, 0.050, 69)
profit_grid = np.array([profit_at_fixed_p(p) for p in p_grid])

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(p_grid, profit_grid, lw=2)
ax.plot(pyo.value(m_a.p), pyo.value(m_a.profit), "o", ms=9, label="bound_frac = 0.01")
ax.plot(pyo.value(m_b.p), pyo.value(m_b.profit), "s", ms=9, label="bound_frac = 0.5")
ax.set_xlabel("Pool composition $p$ [mass fraction]")
ax.set_ylabel("Profit")
ax.set_title("Profit as a function of the pool composition")
ax.legend()
ax.grid(True)
plt.tight_layout()
plt.show()

Comparing the three options¶

summary = pd.DataFrame(
    [
        ["1. Local farms only", "LP", profit_1, "leaves money on the table"],
        ["2. All farms, one truck each", "LP", profit_2, "needs a second truck"],
        ["3. Remote farms pooled", "bilinear NLP", pyo.value(m_a.profit), "nonconvex"],
    ],
    columns=["Option", "Model class", "Profit", "Catch"],
)
display(summary.round(2))

Take away messages¶

Pooling recovers most of the value of the remote farms without a second truck. But the model that says so is nonconvex, and the number it reports depends on where the solver started.
A bilinear term --- a product of two decision variables --- is smooth but irreducibly nonconvex. Unlike $|x|$ , it cannot be reformulated away.
Report the initial point, or the option that determined it. A result you cannot reproduce is a result you cannot defend.
What can be done is to relax the bilinear term: replace it by a new variable bounded by linear functions of its factors. That gives a true upper bound on the profit, and tightening such bounds until they meet the optimum is what the global optimization lecture does.