2  A blend protected against all upper deviations

Problem 02 | Blending with Uncertain Feed Composition | Level 2: Model extension

2.1 Problem Statement

A plant blends 100 t/d of a liquid intermediate. Impurity must not exceed 7 wt%. Costs and available flows are deterministic. Impurity has no beneficial role, so only upward deviations matter for the quality upper bound. There are no volume changes, storage, or blending losses.

Feed Cost (USD/t) Available (t/d) Nominal impurity Maximum upward deviation
A 100 60 0.02 0.01
B 80 60 0.08 0.02
C 55 60 0.14 0.03

All fractions are mass fractions. Deviations of 0.01 mean one percentage point, not a 1% relative change.

Your task. Keep the same production target, prices, and availability. Require the specification to hold for every composition in the box \(a_i\le\tilde a_i\le a_i+d_i\). Find the robust blend and the cost premium relative to Problem 1.

2.1.1 Before you code

Is testing only the average of the upper and lower compositions enough to guarantee quality? Which box corner is worst when all flows are nonnegative?

2.1.2 Deliverables

  1. Define decision variables, units, objective, and constraints. State which data and assumptions are fixed.
  2. Predict one feature of the solution and explain why it should occur.
  3. Build and solve a Pyomo model. Check balances and bounds independently.
  4. Interpret the result and investigate the follow-up question below. For an open-ended challenge, state and defend your chosen policy separately from the reference solution.

2.2 Model Formulation

Let \(x_i\) be feed \(i\) in t/d, with \(0\le x_i\le60\). Define \(c_i\) as cost, \(a_i\) as nominal impurity, and \(d_i\) as maximum upward deviation. Let \(Q=100\) t/d and \(L=0.07\).

The mass balance and cost objective are \[\sum_i x_i=Q,\qquad \min C=\sum_i c_i x_i.\] The impurity constraint compares impurity mass flow with \(LQ\). Multiplying a fixed quality limit by the fixed product rate avoids a ratio.

Because \(x_i\ge0\), the maximum impurity occurs at the upper corner: \[\max_{0\le u_i\le1}\sum_i(a_i+d_i u_i)x_i=\sum_i(a_i+d_i)x_i.\] Thus the robust counterpart is simply \[\sum_i(a_i+d_i)x_i\le LQ.\] No probabilities are needed. The guarantee is conditional on the stated bounds being valid.

2.3 Pyomo Implementation

The code below is the model-building portion of models/blending.py. The level argument selects this problem’s assumptions. Run the complete module from the source bundle to reproduce the solution, including its independent checks:

~/.venvs/optim/bin/python -m models.blending 2
"""Blend design: nominal, box-robust, and budget-robust quality."""

import itertools
import pyomo.environ as pyo
from models.common import value, solve

COST = {"A": 100, "B": 80, "C": 55}  # USD/t
NOMINAL = {"A": 0.02, "B": 0.08, "C": 0.14}  # mass fraction
DELTA = {"A": 0.01, "B": 0.02, "C": 0.03}
CAP = {"A": 60, "B": 60, "C": 60}  # t/d
TOTAL = 100  # t/d
LIMIT = 0.07  # maximum impurity mass fraction


def build(level, gamma=1.5):
    m = pyo.ConcreteModel(name="Blend quality")
    m.I = pyo.Set(initialize=list(COST))
    m.x = pyo.Var(m.I, domain=pyo.NonNegativeReals, bounds=lambda m, i: (0, CAP[i]))
    m.mass = pyo.Constraint(expr=sum(m.x[i] for i in m.I) == TOTAL)
    nominal = sum(NOMINAL[i] * m.x[i] for i in m.I)
    if level == 1:
        m.quality = pyo.Constraint(expr=nominal <= LIMIT * TOTAL)
    elif level == 2:
        m.quality = pyo.Constraint(
            expr=nominal + sum(DELTA[i] * m.x[i] for i in m.I) <= LIMIT * TOTAL
        )
    else:
        m.p = pyo.Var(domain=pyo.NonNegativeReals)
        m.q = pyo.Var(m.I, domain=pyo.NonNegativeReals)
        m.support = pyo.Constraint(
            m.I, rule=lambda m, i: m.p + m.q[i] >= DELTA[i] * m.x[i]
        )
        m.quality = pyo.Constraint(
            expr=nominal + gamma * m.p + sum(m.q[i] for i in m.I) <= LIMIT * TOTAL
        )
    m.cost = pyo.Objective(expr=sum(COST[i] * m.x[i] for i in m.I))
    m._gamma = gamma
    return m

The common solve function calls appsi_highs, checks optimal termination before loading a solution, and checks constraint residuals, bounds, and integer domains. After building the model, use:

from models.common import solve
m = solve(build(2))

2.4 Optimal Solution

The reference objective is 8,771.4286 USD/d. This value is optimal for the explicitly stated model and data.

Feed Flow (t/d)
A 60.0000
B 22.8571
C 17.1429
Quantity Value
Nominal impurity (%) 5.4286
Full-box impurity (%) 7.0000
Protected impurity (%) 7.0000
Gamma 3.0000
Figure 2.1: Blending with Uncertain Feed Composition: reference solution for level 2.

2.4.1 Check your solution

Check all eight corners of the uncertainty box. Calculate the cost premium as \((C_{box}-C_{nominal})/C_{nominal}\), using the nominal optimum as the denominator.

The automated residual, bound, and integrality audit passed with a maximum violation of 0.00e+00 in model units. Domain-specific checks passed. Model units differ across equations, so this numerical audit does not replace dimensional analysis.

2.5 Brief Discussion

The protected design costs 8,771.43 USD/d, a premium of 646.43 USD/d (7.96%) over the nominal optimum. Full-box impurity is 7.000%. Protection at a smaller budget must not be described as full-box protection.

2.5.1 Extension to investigate

Reduce the allowable impurity to 6 wt%. Determine whether the available low-impurity feed still makes the robust model feasible.