19 Chemical Distribution Through Depots

Source and adaptation

Based on Distribution 1, Section 12.19 of Williams (2013). New two-plant, three-depot, four-customer transportation instance. All shipments pass through a depot. The results below solve the stated instance and are not presented as the numerical answer to an unmodified textbook problem.

19.1 Problem Statement

Plants P1 and P2 can supply at most 60 and 70 t/day. Depots D1,D2,D3 have throughput limits 45,40,50 t/day. Customer demands C1–C4 are 20,25,30,20 t/day and must all be met. Plant-to-depot shipping costs (D1,D2,D3) are P1:(2,4,5), P2:(5,2,3) USD/t. Depot-to-customer costs (C1,C2,C3,C4) are D1:(2,3,6,7), D2:(5,2,3,4), D3:(6,5,2,2) USD/t. No losses or inventory are allowed. Minimize total shipping cost.

19.1.1 Inbound costs (USD/t)

Plant	D1	D2	D3	Supply limit (t/day)
P1	2	4	5	60
P2	5	2	3	70

19.1.2 Outbound costs (USD/t)

Depot	C1	C2	C3	C4	Throughput (t/day)
D1	2	3	6	7	45
D2	5	2	3	4	40
D3	6	5	2	2	50

19.1.3 Customer demands

Customer	Demand (t/day)
C1	20
C2	25
C3	30
C4	20

19.2 Model Formulation

Let \(f_{pd}\ge0\) be plant-to-depot flow and \(g_{dc}\ge0\) be depot-to-customer flow in t/day. Minimize \(\sum_{pd}c_{pd}f_{pd}+\sum_{dc}k_{dc}g_{dc}\) subject to

\[\sum_df_{pd}\le U_p,\quad \sum_pf_{pd}=\sum_cg_{dc},\quad\sum_cg_{dc}\le V_d,\quad\sum_dg_{dc}=D_c.\]

The depot balance forbids material creation or disappearance. Plant capacity is an upper bound, while customer demand is an equality. This is a transshipment LP with two transport legs.

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.

19.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 problem19.py. To solve and print this problem independently:

~/.venvs/optim/bin/python -m models.common 19

import pyomo.environ as pyo
from models.common import frame

UP = [60, 70]
VD = [45, 40, 50]
DEMAND = [20, 25, 30, 20]
C = [[2, 4, 5], [5, 2, 3]]
K = [[2, 3, 6, 7], [5, 2, 3, 4], [6, 5, 2, 2]]


def build():
    m = pyo.ConcreteModel()
    m.P = pyo.RangeSet(0, 1)
    m.D = pyo.RangeSet(0, 2)
    m.C = pyo.RangeSet(0, 3)
    m.f = pyo.Var(m.P, m.D, domain=pyo.NonNegativeReals)
    m.g = pyo.Var(m.D, m.C, domain=pyo.NonNegativeReals)
    m.plant = pyo.Constraint(
        m.P, rule=lambda m, p: sum(m.f[p, d] for d in m.D) <= UP[p]
    )
    m.balance = pyo.Constraint(
        m.D,
        rule=lambda m, d: sum(m.f[p, d] for p in m.P) == sum(m.g[d, c] for c in m.C),
    )
    m.cap = pyo.Constraint(m.D, rule=lambda m, d: sum(m.g[d, c] for c in m.C) <= VD[d])
    m.demand = pyo.Constraint(
        m.C, rule=lambda m, c: sum(m.g[d, c] for d in m.D) == DEMAND[c]
    )
    m.obj = pyo.Objective(
        expr=sum(C[p][d] * m.f[p, d] for p in m.P for d in m.D)
        + sum(K[d][c] * m.g[d, c] for d in m.D for c in m.C)
    )
    return m


def check(m):
    assert abs(sum(pyo.value(m.f[p, d]) for p in m.P for d in m.D) - sum(DEMAND)) < 1e-6
    for d in m.D:
        assert (
            abs(
                sum(pyo.value(m.f[p, d]) for p in m.P)
                - sum(pyo.value(m.g[d, c]) for c in m.C)
            )
            < 1e-6
        )


def tables(m):
    return {
        "shipments": frame(
            [
                [f"P{p + 1}", f"D{d + 1}", pyo.value(m.f[p, d])]
                for p in m.P
                for d in m.D
                if pyo.value(m.f[p, d]) > 1e-6
            ]
            + [
                [f"D{d + 1}", f"C{c + 1}", pyo.value(m.g[d, c])]
                for d in m.D
                for c in m.C
                if pyo.value(m.g[d, c]) > 1e-6
            ],
            ["From", "To", "Flow (t/day)"],
        )
    }


def plot(m):
    return (
        ["D1", "D2", "D3"],
        [sum(pyo.value(m.g[d, c]) for c in m.C) for d in m.D],
        "Depot throughput (t/day)",
    )

# 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))

19.4 Optimal Solution

The solver reports optimal termination. The objective is 435 USD/day (minimize). The largest violation across active constraints, variable bounds, and integer domains is 0.00e+00 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.

19.4.1 Shipments

From	To	Flow (t/day)
P1	D1	25
P2	D2	40
P2	D3	30
D1	C1	20
D1	C2	5
D2	C2	20
D2	C3	20
D3	C3	10
D3	C4	20

Figure 19.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.

19.5 Brief Discussion

The network delivers 95 t/day at a total transport cost of 435.00 USD/day, or 4.58 USD/t averaged over deliveries.

A cheap last-mile route can be unattractive after including inbound depot cost. Depot capacity can redirect flows across the whole network. Experiment: add 10 t/day capacity at each depot separately and report which expansion reduces cost the most.