13  Balanced Assignment of Chemical Distributors

NoteSource and adaptation

Based on Market sharing, Section 12.13 of Williams (2013). New eight-distributor instance with three balancing measures and a 15 percentage-point feasibility tolerance. The results below solve the stated instance and are not presented as the numerical answer to an unmodified textbook problem.

13.1 Problem Statement

Allocate each of eight distributors entirely to division D1 or D2. Distributor records (solvent volume, resin volume, delivery points) are (12,8,3), (8,15,4), (15,6,5), (10,12,3), (6,14,2), (14,9,4), (9,11,3), (16,5,6). Volumes are thousand t/year and points are counts. D1 should receive 40% of each total; D2 receives the remainder. Every D1 share must lie between 25% and 55%. Minimize the maximum absolute deviation from the 40% target across the three measures. A distributor cannot be split.

13.1.1 Distributor data

Distributor Solvent (1000 t/year) Resin (1000 t/year) Delivery points
1 12 8 3
2 8 15 4
3 15 6 5
4 10 12 3
5 6 14 2
6 14 9 4
7 9 11 3
8 16 5 6

13.2 Model Formulation

Let \(x_i\in\{0,1\}\) indicate assignment to D1. For measure \(k\), known normalized weights are \(w_{ik}=a_{ik}/\sum_j a_{jk}\). Let \(E\ge0\).

\[0.25\le\sum_i w_{ik}x_i\le0.55,\qquad -E\le\sum_iw_{ik}x_i-0.4\le E.\]

Minimize \(E\). Normalizing each criterion prevents physical units or large totals from dominating the comparison. The result is a minimax MILP; \(E=0.01\) means one percentage point, not one percent relative to the target.

All decision variables and units refer to the problem statement above. Continuous variables are nonnegative unless explicitly stated otherwise; binary and integer domains are specified in the equations and code.

13.3 Pyomo Implementation

The following Python implementation uses Pyomo and HiGHS. Run from the workbook root so that the models package and shared helpers are importable. The shared solver and audit functions check optimal termination before loading values, then verify every active constraint, variable bound, and integer domain. Any imported earlier-chapter model supplies the data and balances already explained there.

The full source is problem13.py. To solve and print this problem independently:

~/.venvs/optim/bin/python -m models.common 13
import itertools
import pyomo.environ as pyo
from models.common import frame

DATA = [
    (12, 8, 3),
    (8, 15, 4),
    (15, 6, 5),
    (10, 12, 3),
    (6, 14, 2),
    (14, 9, 4),
    (9, 11, 3),
    (16, 5, 6),
]
TOTAL = [sum(row[k] for row in DATA) for k in range(3)]


def build():
    m = pyo.ConcreteModel()
    m.I = pyo.RangeSet(0, 7)
    m.x = pyo.Var(m.I, domain=pyo.Binary)
    m.E = pyo.Var(domain=pyo.NonNegativeReals)
    m.c = pyo.ConstraintList()
    for k in range(3):
        share = sum(DATA[i][k] / TOTAL[k] * m.x[i] for i in m.I)
        m.c.add(share >= 0.25)
        m.c.add(share <= 0.55)
        m.c.add(share - 0.4 <= m.E)
        m.c.add(0.4 - share <= m.E)
    m.obj = pyo.Objective(expr=m.E)
    return m


def check(m):
    candidates = []
    for bits in itertools.product([0, 1], repeat=8):
        shares = [
            sum(DATA[i][k] * bits[i] for i in range(8)) / TOTAL[k] for k in range(3)
        ]
        if all(0.25 <= s <= 0.55 for s in shares):
            candidates.append(max(abs(s - 0.4) for s in shares))
    assert abs(pyo.value(m.obj) - min(candidates)) < 1e-6


def tables(m):
    return {
        "assignment": frame(
            [[i + 1, "D1" if pyo.value(m.x[i]) > 0.5 else "D2"] for i in m.I],
            ["Distributor", "Division"],
        ),
        "shares": frame(
            [
                [
                    name,
                    100 * sum(DATA[i][k] * pyo.value(m.x[i]) for i in m.I) / TOTAL[k],
                ]
                for k, name in enumerate(["Solvent", "Resin", "Delivery points"])
            ],
            ["Measure", "D1 share (%)"],
        ),
    }


def plot(m):
    return (
        ["Solvent", "Resin", "Delivery points"],
        [
            100 * sum(DATA[i][k] * pyo.value(m.x[i]) for i in m.I) / TOTAL[k]
            for k in range(3)
        ],
        "D1 share (%)",
    )

# Solve, audit constraints, and run domain-specific checks.
from models.common import solve, audit
model = solve(build())
audit(model)
check(model)
for name, result_table in tables(model).items():
    print(name)
    print(result_table.to_string(index=False))

13.4 Optimal Solution

The solver reports optimal termination. The objective is 0.033333 maximum absolute share deviation (fraction) (minimize). The largest violation across active constraints, variable bounds, and integer domains is 2.78e-17 in the corresponding model units. The problem-specific checks also pass. These checks establish numerical consistency with the stated model, not the validity of its assumptions for a real facility.

13.4.1 Assignment

Distributor Division
1 D2
2 D1
3 D2
4 D2
5 D2
6 D2
7 D1
8 D1

13.4.2 Shares

Measure D1 share (%)
Solvent 36.6667
Resin 38.75
Delivery points 43.3333
Figure 13.1: Selected quantities from the verified optimal solution. Units are stated on the axis.

Tables round numerical values for reading; feasibility checks use the original solver values. Multiple optimal decisions may exist. Machine-readable result records the solver status and package versions. Figure-generation code is in figures/workbook.py.

13.5 Brief Discussion

The best worst-case deviation is 3.333 percentage points from the 40% target. This is a deviation in share, not in physical shipment quantity.

Indivisibility usually prevents every measure from attaining exactly 40%. A small deviation in one metric does not prove the allocation is balanced in the others. Experiment: minimize the sum of absolute deviations and compare which criterion becomes worst.