Product allocation and functional-area sizing
Aims
Determine the allocation of products to material flows and the sizes of the cross-docking, reserve, and forward areas.
Learning objectives
After completing this tutorial, you will be able to:
- calculate order quantities and average dwell times;
- formulate a product-flow allocation and functional-area sizing model; and
- interpret product assignments, annual costs, and space utilization.
Product allocation and functional-area sizing
A company is choosing the allocation of products and the sizes of functional areas to minimize total annual handling and storage costs. The warehouse manages six product classes, each with unique demand patterns, handling requirements, and storage characteristics. Use the allocation model of Heragu et al. (2005) with unit-load simplifications.
Product flows
When products are assigned to flows 1 (cross-docking), 2 (reserve), or 4 (forward), they are handled as unit loads such as pallets. When assigned to flow 3 (reserve to forward), products are initially stored in the reserve area as unit loads and are later broken down into cartons. Break-bulk adds handling and can require more storage positions. The break-bulk factor \(q_{ij}\) adjusts handling and storage costs for product \(i\) in flow \(j\). For simplicity, we assume \(q_{ij}=1\) for all products and flows in this case study. The worked example therefore uses a common unit-load basis. The formulation below denotes flows by \(f\). Table 18.1 shows the material paths and the areas they use.
| Flow | Material path | Areas using modeled storage capacity |
|---|---|---|
| CD (1) | Receiving → Cross-docking → Shipping | Cross-docking |
| R (2) | Receiving → Reserve → Shipping | Reserve |
| RF (3) | Receiving → Reserve → Forward → Shipping | Reserve and forward |
| F (4) | Receiving → Forward → Shipping | Forward |
Receiving and shipping are boundary interfaces in this model. The RF path includes reserve-to-forward movement under these assumptions.
Product and cost inputs
Table 18.2 shows the costs and space requirements for each product. Annual demand, order cost, price, and floor space per unit load are denoted by \(D_i\), \(K_i\), \(P_i\), and \(s_i\). Only products 3 and 6 have a reserve dwell fraction greater than zero. They spend on average 20% and 100% of their time in the reserve area, respectively, when they are assigned to flow 3 (reserve to forward). The reserve dwell share \(\rho_i\) applies only if the product uses RF.
| Product | Loads/year | Order cost ($) | Price ($/load) | m²/load | Reserve share in RF |
|---|---|---|---|---|---|
| 1 | 10,000 | 50 | 500 | 10 | 0.00 |
| 2 | 15,000 | 50 | 650 | 15 | 0.00 |
| 3 | 25,000 | 50 | 350 | 25 | 0.20 |
| 4 | 2,000 | 50 | 250 | 10 | 0.00 |
| 5 | 1,500 | 50 | 225 | 12 | 0.00 |
| 6 | 95,000 | 150 | 150 | 13 | 1.00 |
Table 18.3 shows the handling cost \(h_{if}\) per unit load per transaction for each product in each material flow. Table 18.4 shows the yearly cost \(c_{if}\) of storing an average product unit load in each material flow. Both costs depend on the assigned flow \(f\).
| Product | CD | R | RF | F |
|---|---|---|---|---|
| 1 | 0.0707 | 0.0849 | 0.1061 | 0.0778 |
| 2 | 0.0203 | 0.2023 | 0.2023 | 0.2023 |
| 3 | 0.0267 | 0.0420 | 0.0054 | 0.0481 |
| 4 | 0.3354 | 0.5590 | 1.0062 | 0.0671 |
| 5 | 0.4083 | 0.6804 | 1.2248 | 0.8165 |
| 6 | 0.0726 | 0.0871 | 0.1088 | 0.0798 |
| Product | CD | R | RF | F |
|---|---|---|---|---|
| 1 | 20 | 5 | 10 | 15 |
| 2 | 15 | 5 | 10 | 10 |
| 3 | 4 | 20 | 1 | 9 |
| 4 | 5 | 4 | 5 | 1 |
| 5 | 15 | 25 | 45 | 30 |
| 6 | 20 | 5 | 10 | 15 |
The warehouse has 100,000 m² of storage space available, which must be allocated among cross-docking, reserve, and forward storage. Table 18.5 gives the lower and upper bounds for each functional area. The lower bounds can require substantially more floor area than the modeled average inventory occupies. Each functional area has 1 vertical storage level.
| Area | Minimum (m²) | Maximum (m²) | Levels |
|---|---|---|---|
| CD | 0 | 15,000 | 1 |
| R | 35,000 | 75,000 | 1 |
| F | 35,000 | 75,000 | 1 |
We load the product, handling-cost, storage-cost, and area-bound tables from CSV files in tutorial-inputs/. The function below also calculates EOQ, average inventory, and average dwell time. Its complete source and the CSV contents are shown under Run the Python example at the end of this tutorial.
EOQ and average dwell-time calculation
We calculate the Economic Order Quantity (EOQ) and the average inventory level for each product, using the EOQ formula. With annual carrying rate \(r=0.10\), the holding cost used in EOQ is \(H_i=rP_i\). Under the example’s steady demand and instantaneous receipt assumptions:
\[ Q_i=\sqrt{\frac{2D_iK_i}{H_i}},\qquad \bar I_i=\frac{Q_i}{2},\qquad T_i=\frac{\bar I_i}{D_i}. \tag{18.1}\]
\(Q_i\) is the purchasing order quantity, \(\bar I_i\) is average inventory in unit loads, and \(T_i\) is average dwell time in years. The time between two successive purchasing receipts is \(Q_i/D_i\). For example, if the order quantity is 100 unit loads and the annual demand is 1,000 unit loads, the time between receipts is 0.1 years, or approximately 36.5 days. The average inventory level is half the order quantity, which in this case is 50 unit loads. The average dwell time is \(100/(2\times1{,}000)=0.05\) years, or approximately 18.25 days.
Table 18.6 shows the calculated values for the 6 products.
| Product | EOQ (loads) | Average inventory (loads) | Dwell time (years) |
|---|---|---|---|
| 1 | 141.42 | 70.71 | 0.007071 |
| 2 | 151.91 | 75.96 | 0.005064 |
| 3 | 267.26 | 133.63 | 0.005345 |
| 4 | 89.44 | 44.72 | 0.022361 |
| 5 | 81.65 | 40.82 | 0.027217 |
| 6 | 1,378.40 | 689.20 | 0.007255 |
The carrying cost \(H_i\) sets the EOQ. The flow-dependent storage cost \(c_{if}\) enters the allocation objective. Keep those two costs distinct. \(Q_i\) and \(\rho_i\) are calculated or supplied before allocation. The allocation model does not choose forward replenishment quantities or frequency. Average inventory satisfies \(Q_i/2=D_iT_i\), so either expression gives the inventory term in the objective (Heragu et al. 2005).
Calculate occupied space using average inventory \(Q_i/2\). Average space does not establish that peak inventory, purchasing receipts, or safety stock will fit. Check those conditions separately before implementing a design. The synthetic CD flow also uses the supplied average-inventory approximation; a real cross-dock needs evidence about staging dwell and peak arrivals.
Model formulation
The goal is to determine the allocation of products to material flows and the sizing of functional areas while adhering to constraints on space availability and product-specific requirements. We determine the product assignments and the sizes of the cross-docking, reserve, and forward areas together.
Let \(x_{if}=1\) when product \(i\) uses flow \(f\), and zero otherwise. Let \(a_{\text{CD}},a_{\text{R}},a_{\text{F}}\) be allocated floor areas in m². An equivalent formulation uses area fractions \(w_a\), with \(a_a=100{,}000w_a\) for this one-level case.
The objective function minimizes total annual handling and storage costs. The baseline convention counts two handling transactions for every flow:
\[ \min \sum_i\sum_{f=1}^{4} \left(2D_i h_{if}+\bar I_i c_{if}\right)x_{if}. \tag{18.2}\]
The objective is annual handling plus storage cost in dollars/year. Area allocations enter the capacity constraints; there is no extra floor-area cost term.
Each product must be assigned to exactly one flow:
\[ \sum_{f=1}^{4}x_{if}=1\quad(\forall i),\qquad x_{if}\in\{0,1\}. \]
The total space used in each functional area must not exceed the allocated space:
\[ \begin{aligned} \sum_i \bar I_i s_i x_{i1} &\le a_{\text{CD}},\\ \sum_i \bar I_i s_i(x_{i2}+\rho_i x_{i3}) &\le a_{\text{R}},\\ \sum_i \bar I_i s_i((1-\rho_i)x_{i3}+x_{i4}) &\le a_{\text{F}}. \end{aligned} \tag{18.3}\]
The total allocated space must equal the total available space. The size of each functional area must be within its specified lower and upper bounds:
\[ a_{\text{CD}}+a_{\text{R}}+a_{\text{F}}=100{,}000,\qquad L_a\le a_a\le U_a\quad(a\in\{\text{CD},\text{R},\text{F}\}). \]
For RF, changing the product-flow assignment changes occupancy in both reserve and forward. The supplied \(\rho_i\) divides that product’s average inventory between them. The cost coefficients and dwell shares remain fixed when the forward allocation changes. The model therefore cannot calculate a reduction in replenishment work from a larger forward allocation alone.
Model implementation
First, we create the optimization model and define the decision variables. Then we add the objective and the constraints in PuLP:
Show code
import pulp as pl
model = pl.LpProblem("AllocationSizing", pl.LpMinimize)
x = pl.LpVariable.dicts("x", costs, cat=pl.LpBinary)
area = {
a: pl.LpVariable(f"area_{a}", bounds.loc[a, "Lower"],
bounds.loc[a, "Upper"])
for a in bounds.index
}
model += pl.lpSum(costs[key] * x[key] for key in costs)
for i in products.index:
model += pl.lpSum(x[i, f] for f in (1, 2, 3, 4)) == 1
for a in bounds.index:
model += pl.lpSum(
occupied[key][a] * x[key] for key in costs
) <= area[a]
model += pl.lpSum(area.values()) == 100000The Python script solves this formulation and saves the results. With only 6 products and 4 flows, you can check the optimum by enumerating all \(4^6=4{,}096\) assignments.
Product allocation and area sizes
Table 18.7 reports the chosen flows and cost components. Table 18.8 reports one corresponding feasible area allocation.
| Product | Flow | Handling ($/year) | Storage ($/year) | Total ($/year) |
|---|---|---|---|---|
| 1 | R | 1,698.00 | 353.55 | 2,051.55 |
| 2 | CD | 609.00 | 1,139.33 | 1,748.33 |
| 3 | RF | 270.00 | 133.63 | 403.63 |
| 4 | F | 268.40 | 44.72 | 313.12 |
| 5 | CD | 1,224.90 | 612.37 | 1,837.27 |
| 6 | R | 16,549.00 | 3,446.01 | 19,995.01 |
| Total | All | 20,619.30 | 5,729.62 | 26,348.92 |
| Area | Occupied (m²) | Allocated (m²) | Utilization (%) |
|---|---|---|---|
| CD | 1,629.23 | 15,000.00 | 10.86 |
| R | 10,334.89 | 35,000.00 | 29.53 |
| F | 3,119.83 | 50,000.00 | 6.24 |
The allocated space sums to 100,000 m². The occupied space is much smaller because the case’s large area lower bounds and full-allocation equality leave considerable slack. Several area splits achieve the same minimum cost. For example, transferring 15,000 m² from forward to reserve gives \((a_{\text{CD}},a_{\text{R}},a_{\text{F}})=(15{,}000,50{,}000,35{,}000)\) with the same assignments and cost. Check both capacity and bounds to confirm this statement.
Alternative handling convention
The baseline counts 2 handling transactions for every flow. An alternative counts 3 for RF, covering receipt, reserve-to-forward transfer, and shipping. To evaluate that alternative, run:
The choice must match the definition of the handling coefficient. If that coefficient already includes the whole path, multiplying it again would double-count handling. Record the convention and compare the resulting assignments and cost components.
Using warehouse data
Reserve full-case picks and forward picks may serve the same SKU. Justify any product-class definition or exclusive-flow simplification before using \(\sum_f x_{if}=1\). Purchasing costs, dwell shares, and cross-docking eligibility need supporting data. Document each missing input and separate assumed scenarios from company facts.
RF uses space in both reserve and forward, and its handling includes the transfer between them. The allocation model calculates that space and cost using fixed dwell shares and flow-specific costs. Check replenishment access when evaluating whether the area allocations are feasible. Additional claims about changing replenishment frequency require supporting calculations and evidence.
Save all numerical outputs and the complete model settings with the code.
Use Python 3.11 or later. Save the source and input files below using the paths shown above each code block. Run the snippets in order from the directory containing scripts/ and tutorial-inputs/. The allocation model uses PuLP with the CBC solver.
The first command installs the required packages. The second solves the model and writes its results to tutorial-results/. The third checks the saved assignments, costs, and area capacities without solving again.
Show complete Python source
scripts/layout_tutorial.py
"""Reproduce the worked storage allocation and warehouse layout examples.
Install pandas, matplotlib, tabulate, pulp, gurobipy and the pinned whlayout
package listed in ../tutorial-requirements.txt. --verify needs only pandas.
"""
from __future__ import annotations
import argparse
import hashlib
import itertools
import json
import math
import time
from pathlib import Path
import pandas as pd
ROOT = Path(__file__).resolve().parents[1]
AREA_NAMES = ("CD", "R", "F")
FLOW_NAMES = {1: "CD", 2: "R", 3: "RF", 4: "F"}
WHLAYOUT_REVISION = "c2e789ce76db9aafae385f6ab2e5ce78b652cc6d"
def source_digest(path):
"""Hash text with LF newlines so Windows and Linux checkouts agree."""
return hashlib.sha256(path.read_bytes().replace(b"\r\n", b"\n")).hexdigest()
def load_allocation_inputs():
"""Read the six-product historical case, with stable product/flow IDs."""
source = ROOT / "tutorial-inputs"
products = pd.read_csv(source / "products.csv", index_col="Product")
handling = pd.read_csv(source / "handling.csv", index_col="Product")
storage = pd.read_csv(source / "storage.csv", index_col="Product")
bounds = pd.read_csv(source / "area-bounds.csv", index_col="Area")
assert products.index.equals(handling.index)
assert products.index.equals(storage.index)
products["EOQ"] = (
2 * products.Demand * products.OrderCost
/ (products.Price * products.CarryingRate)
) ** 0.5
products["AverageInventory"] = products.EOQ / 2
products["DwellYears"] = products.AverageInventory / products.Demand
return products, handling, storage, bounds
def allocation_coefficients(products, handling, storage, rf_operations=2):
"""Historical q=1, one level, average-inventory allocation model."""
costs, occupancy = {}, {}
for i, row in products.iterrows():
for f in FLOW_NAMES:
operations = rf_operations if f == 3 else 2
costs[i, f] = (
operations * row.Demand * handling.loc[i, str(f)]
+ row.AverageInventory * storage.loc[i, str(f)]
)
base = row.AverageInventory * row.SpacePerLoad
occupancy[i, f] = {
"CD": base if f == 1 else 0.0,
"R": base * (1 if f == 2 else row.ReserveShare
if f == 3 else 0),
"F": base * (1 if f == 4 else 1 - row.ReserveShare
if f == 3 else 0),
}
return costs, occupancy
def validate_allocation(assignments, areas, bounds, total=100000.0):
"""Verify independent capacity and accounting conditions."""
assert len(assignments) == 6
assert assignments.Product.nunique() == 6
assert assignments.FlowID.isin(FLOW_NAMES).all()
assert abs(areas.Allocated.sum() - total) < 1e-5
for row in areas.itertuples():
assert bounds.loc[row.Area, "Lower"] - 1e-5 <= row.Allocated
assert row.Allocated <= bounds.loc[row.Area, "Upper"] + 1e-5
assert row.Occupied <= row.Allocated + 1e-5
assert row.Occupied >= 0
def solve_allocation(output, rf_operations=2):
import pulp as pl
products, handling, storage, bounds = load_allocation_inputs()
costs, occupancy = allocation_coefficients(
products, handling, storage, rf_operations
)
model = pl.LpProblem("HistoricalAllocationSizing", pl.LpMinimize)
x = pl.LpVariable.dicts("x", costs, cat=pl.LpBinary)
area = {
a: pl.LpVariable(f"area_{a}", bounds.loc[a, "Lower"],
bounds.loc[a, "Upper"])
for a in AREA_NAMES
}
model += pl.lpSum(costs[key] * x[key] for key in costs)
for i in products.index:
model += pl.lpSum(x[i, f] for f in FLOW_NAMES) == 1
used = {
a: pl.lpSum(occupancy[key][a] * x[key] for key in costs)
for a in AREA_NAMES
}
for a in AREA_NAMES:
model += used[a] <= area[a]
model += pl.lpSum(area.values()) == 100000.0
started = time.perf_counter()
model.solve(pl.PULP_CBC_CMD(msg=False, threads=1, timeLimit=30))
status = pl.LpStatus[model.status]
metadata = {
"status": status, "runtime_seconds": time.perf_counter() - started,
"solver": "PuLP/CBC", "rf_handling_operations": rf_operations,
"total_area_m2": 100000.0, "break_bulk_factor": 1,
"vertical_levels": 1, "inventory_basis": "EOQ/2 average inventory",
}
(output / "allocation-run.json").write_text(
json.dumps(metadata, indent=2), encoding="utf-8"
)
if status != "Optimal":
raise RuntimeError(f"Allocation solve returned {status}")
rows = []
for i in products.index:
f = next(f for f in FLOW_NAMES if pl.value(x[i, f]) > 0.5)
operations = rf_operations if f == 3 else 2
rows.append({
"Product": i, "FlowID": f, "Flow": FLOW_NAMES[f],
"HandlingCost": operations * products.loc[i, "Demand"]
* handling.loc[i, str(f)],
"StorageCost": products.loc[i, "AverageInventory"]
* storage.loc[i, str(f)],
})
assignments = pd.DataFrame(rows)
assignments["TotalCost"] = (
assignments.HandlingCost + assignments.StorageCost
)
areas = pd.DataFrame([
{"Area": a, "Occupied": pl.value(used[a]),
"Allocated": pl.value(area[a]),
"UtilizationPct": 100 * pl.value(used[a]) / pl.value(area[a])
if pl.value(area[a]) > 0 else 0.0}
for a in AREA_NAMES
])
validate_allocation(assignments, areas, bounds)
assert math.isclose(
assignments.TotalCost.sum(), pl.value(model.objective), abs_tol=1e-5
)
# Independent enumeration of the 4^6 assignments validates this tiny MILP.
valid_costs = []
for flows in itertools.product(FLOW_NAMES, repeat=len(products)):
allocation = dict(zip(products.index, flows))
occupied = {
a: sum(occupancy[i, f][a] for i, f in allocation.items())
for a in AREA_NAMES
}
if all(occupied[a] <= bounds.loc[a, "Upper"] for a in AREA_NAMES):
required = sum(max(occupied[a], bounds.loc[a, "Lower"])
for a in AREA_NAMES)
if required <= 100000 <= bounds.Upper.sum():
valid_costs.append(sum(costs[i, f]
for i, f in allocation.items()))
assert math.isclose(min(valid_costs), pl.value(model.objective),
abs_tol=1e-5)
products.to_csv(output / "inventory.csv")
assignments.to_csv(output / "assignments.csv", index=False)
areas.to_csv(output / "areas.csv", index=False)
metadata.update(objective=pl.value(model.objective),
enumeration_check=True, checks_passed=True)
(output / "allocation-run.json").write_text(
json.dumps(metadata, indent=2), encoding="utf-8"
)
def load_layout_inputs():
"""Use the original case legend and its identical staging relationships."""
source = ROOT / "instances/data/coolblue"
raw = pd.read_csv(source / "adjacency_matrix.csv", index_col=0)
legend = pd.read_csv(source / "muther_legend.csv", index_col=0)
areas = pd.read_csv(source / "space_program_abs.csv", index_col=0).Sq_m
staging = ["Outbound Staging - 2-Man Delivery",
"Outbound Staging - Parcel"]
selected = ["Inbound Dock", "Receiving/Staging",
"Pallet Reserve Storage (Bulk)", "Packing / Wrap / Banding",
"Shipping Dock"]
assert raw.loc[staging[0], selected].equals(
raw.loc[staging[1], selected]
), "The historical aggregation requires identical staging relationships."
weights = raw.replace("-", "U").apply(
lambda column: column.map(legend.Weight.to_dict())
).astype(float)
weights = weights.loc[selected + [staging[0]], selected + [staging[0]]]
weights = weights.rename(index={staging[0]: "Outbound Staging"},
columns={staging[0]: "Outbound Staging"})
areas = areas.loc[selected].to_dict() | {
"Outbound Staging": float(areas.loc[staging].sum())
}
assert (weights.values >= 0).all()
assert (weights.values == weights.values.T).all()
for d in weights.index:
weights.loc[d, d] = 0
return areas, weights
def layout_configurations(areas, weights):
from whlayout.layoututils import propose_building_sides
bx, by = propose_building_sides(areas, slack=1.0, aspect_ratio=1.5)
configurations = []
for shape in ("I", "L", "U"):
current_areas = dict(areas)
current_weights = weights.copy()
if shape == "L":
current_areas["L-Corner"] = 0.25 * bx * by
current_weights["L-Corner"] = 0.0
current_weights.loc["L-Corner"] = 0.0
size = propose_building_sides(
current_areas, slack=1.0, aspect_ratio=1.5
)
anchors = {
"Inbound Dock": [{"side": "south", "offset": 0}],
"Shipping Dock": [{"side": "east", "offset": 0}],
"L-Corner": [{"side": "east", "offset": 0},
{"side": "south", "offset": 0}],
}
elif shape == "I":
size = (bx, by)
anchors = {
"Inbound Dock": [{"side": "west", "offset": 0}],
"Shipping Dock": [{"side": "east", "offset": 0}],
}
else:
size = (bx, by)
anchors = {
"Inbound Dock": [{"side": "south", "offset": 0}],
"Shipping Dock": [{"side": "south", "offset": 0}],
}
skeleton = []
for d in areas:
if d in anchors:
continue
if shape == "I":
skeleton.extend([(d, "east", "Inbound Dock"),
(d, "west", "Shipping Dock")])
elif shape == "L":
skeleton.extend([(d, "north", "Inbound Dock"),
(d, "west", "Shipping Dock")])
else:
skeleton.extend([(d, "north", "Inbound Dock"),
(d, "north", "Shipping Dock")])
configurations.append({
"shape": shape, "areas": current_areas,
"weights": current_weights, "building": size,
"anchors": anchors, "skeleton": skeleton,
})
return configurations
def validate_geometry(table, bx, by, anchors, skeleton, tol=1e-3):
"""Validate exported geometry independently of solver helper variables."""
assert table.index.is_unique
assert table[["width", "height"]].gt(0).all().all()
assert (table.x_west >= -tol).all()
assert (table.y_south >= -tol).all()
assert (table.x_east <= bx + tol).all()
assert (table.y_north <= by + tol).all()
assert ((table.width - (table.x_east - table.x_west)).abs() <= tol).all()
assert ((table.height - (table.y_north - table.y_south)).abs() <= tol).all()
assert ((table.center_x - (table.x_east + table.x_west) / 2)
.abs() <= tol).all()
assert ((table.center_y - (table.y_north + table.y_south) / 2)
.abs() <= tol).all()
area = table.width * table.height
assert ((area - table.area_target).abs()
<= tol + 1e-6 * table.area_target).all()
for a, b in itertools.combinations(table.index, 2):
p, q = table.loc[a], table.loc[b]
assert (p.x_east <= q.x_west + tol
or q.x_east <= p.x_west + tol
or p.y_north <= q.y_south + tol
or q.y_north <= p.y_south + tol)
for d, specifications in anchors.items():
for specification in specifications:
side = specification["side"]
coordinate, target = {
"west": ("x_west", 0), "east": ("x_east", bx),
"south": ("y_south", 0), "north": ("y_north", by),
}[side]
assert abs(table.loc[d, coordinate] - target) <= tol
for a, direction, b in skeleton:
p, q = table.loc[a], table.loc[b]
assert {
"east": p.x_west >= q.x_east - tol,
"west": p.x_east <= q.x_west + tol,
"north": p.y_south >= q.y_north - tol,
"south": p.y_north <= q.y_south + tol,
}[direction]
def centroid_score(table, weights):
return sum(
weights.loc[a, b] * (
abs(table.loc[a, "center_x"] - table.loc[b, "center_x"])
+ abs(table.loc[a, "center_y"] - table.loc[b, "center_y"])
) for a, b in itertools.combinations(table.index, 2)
)
def draw_layout(table, bx, by, shape, output):
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import textwrap
fig, ax = plt.subplots(figsize=(10, 8), layout="constrained")
colors = plt.get_cmap("Pastel1")
labels = []
for index, (department, row) in enumerate(table.iterrows()):
corner = department == "L-Corner"
ax.add_patch(Rectangle(
(row.x_west, row.y_south), row.width, row.height,
facecolor="#d9dde3" if corner else colors(index % 9),
edgecolor="#24344b", linewidth=1.2,
hatch="//" if corner else None,
))
label = "Void" if corner else f"D{index + 1}"
ax.text(row.center_x, row.center_y, label, ha="center",
va="center", fontsize=12, weight="bold")
labels.append(f"{label}: {department} ({row.width:.1f} × "
f"{row.height:.1f} m; {row.area_target:,.0f} m²)")
ax.set(xlim=(0, bx), ylim=(0, by), xlabel="x (m)", ylabel="y (m)",
title=f"{shape}-shaped flow pattern: Coolblue case")
ax.set_aspect("equal")
caption = "\n".join(textwrap.fill(line, 90) for line in labels)
fig.text(0.08, -0.01, caption, fontsize=10, va="top")
fig.savefig(output / f"layout-{shape.lower()}.svg", bbox_inches="tight")
plt.close(fig)
def solve_layouts(output, time_limit, seed):
# Reuse the historical package, including its exact-area formulation.
from whlayout.model import build_block_layout_model
from whlayout.utils import extract_solution_table
areas, weights = load_layout_inputs()
pd.Series(areas, name="Area_m2").rename_axis("Department").to_csv(
output / "departments.csv"
)
weights.to_csv(output / "weights.csv")
runs = []
for config in layout_configurations(areas, weights):
shape = config["shape"]
bx, by = config["building"]
model, variables = build_block_layout_model(
dept_list=list(config["areas"]), areas_m2=config["areas"],
weights=config["weights"], building_x=bx, building_y=by,
anchors=config["anchors"], design_skeleton=config["skeleton"],
area_calculation="exact", mip_gap=0.0001,
time_limit=time_limit, log_to_console=False,
write_lp=str(output / f"layout-{shape.lower()}.lp"),
)
model.Params.Threads = 1
model.Params.Seed = seed
model.optimize()
run = {
"Shape": shape, "Status": int(model.Status),
"SolutionCount": int(model.SolCount), "Runtime": model.Runtime,
"BuildingX": bx, "BuildingY": by,
"GrossEnvelope": bx * by,
"ProductiveArea": sum(areas.values()),
"CornerVoid": config["areas"].get("L-Corner", 0),
"TimeLimit": time_limit, "Seed": seed, "Threads": 1,
"RequestedGap": 0.0001, "AreaMethod": "exact",
"WhlayoutRevision": WHLAYOUT_REVISION,
"Anchors": config["anchors"], "Skeleton": config["skeleton"],
}
runs.append(run)
(output / "layout-runs.json").write_text(
json.dumps(runs, indent=2), encoding="utf-8"
)
if model.SolCount == 0:
model.dispose()
raise RuntimeError(f"{shape}: no feasible incumbent; run retained")
table = extract_solution_table(
list(config["areas"]), bx, by, variables,
areas_m2=config["areas"],
)
validate_geometry(table, bx, by, config["anchors"], config["skeleton"])
score = centroid_score(table, config["weights"])
assert math.isclose(score, model.ObjVal, rel_tol=1e-6, abs_tol=1e-3)
run.update(Objective=score, Bound=model.ObjBound, Gap=model.MIPGap,
MaxAreaErrorPct=table.area_gap_pct.abs().max(),
ChecksPassed=True)
table.to_csv(output / f"layout-{shape.lower()}.csv")
draw_layout(table, bx, by, shape, output)
(output / "layout-runs.json").write_text(
json.dumps(runs, indent=2), encoding="utf-8"
)
model.dispose()
pd.DataFrame(runs).to_csv(output / "layout-comparison.csv", index=False)
def verify_saved(output, require_layouts=True):
products, handling, storage, bounds = load_allocation_inputs()
assignments = pd.read_csv(output / "assignments.csv")
areas = pd.read_csv(output / "areas.csv")
run = json.loads((output / "allocation-run.json").read_text())
costs, occupancy = allocation_coefficients(
products, handling, storage, run["rf_handling_operations"]
)
validate_allocation(assignments, areas, bounds)
assert run["status"] == "Optimal" and run["checks_passed"]
selected = dict(zip(assignments.Product, assignments.FlowID))
for row in assignments.itertuples():
assert row.Flow == FLOW_NAMES[row.FlowID]
operations = run["rf_handling_operations"] if row.FlowID == 3 else 2
expected_handling = (operations * products.loc[row.Product, "Demand"]
* handling.loc[row.Product, str(row.FlowID)])
expected_storage = (products.loc[row.Product, "AverageInventory"]
* storage.loc[row.Product, str(row.FlowID)])
assert math.isclose(row.HandlingCost, expected_handling)
assert math.isclose(row.StorageCost, expected_storage)
assert math.isclose(row.TotalCost, expected_handling + expected_storage)
assert math.isclose(sum(costs[key] for key in selected.items()),
run["objective"], abs_tol=1e-5)
for row in areas.itertuples():
assert math.isclose(
sum(occupancy[key][row.Area] for key in selected.items()),
row.Occupied, abs_tol=1e-5
)
expected_utilization = (100 * row.Occupied / row.Allocated
if row.Allocated else 0)
assert math.isclose(row.UtilizationPct, expected_utilization)
inventory = pd.read_csv(output / "inventory.csv", index_col="Product")
for column in ("EOQ", "AverageInventory", "DwellYears"):
assert (abs(inventory[column] - products[column]) < 1e-8).all()
layout_file = output / "layout-runs.json"
if require_layouts:
assert layout_file.exists(), "The published tutorial needs all layouts."
if layout_file.exists():
runs = json.loads(layout_file.read_text())
assert {run["Shape"] for run in runs} == {"I", "L", "U"}
areas, weights = load_layout_inputs()
comparison = pd.read_csv(output / "layout-comparison.csv")
for run in runs:
shape = run["Shape"]
current_weights = weights.copy()
current_areas = dict(areas)
if shape == "L":
current_weights["L-Corner"] = 0.0
current_weights.loc["L-Corner"] = 0.0
current_areas["L-Corner"] = run["CornerVoid"]
table = pd.read_csv(output / f"layout-{shape.lower()}.csv",
index_col="Department")
assert set(table.index) == set(current_areas)
for d, target in current_areas.items():
assert math.isclose(table.loc[d, "area_target"], target)
validate_geometry(table, run["BuildingX"], run["BuildingY"],
run["Anchors"], run["Skeleton"])
assert math.isclose(centroid_score(table, current_weights),
run["Objective"], abs_tol=1e-3)
reported = comparison.loc[comparison.Shape == shape].iloc[0]
assert math.isclose(reported.Objective, run["Objective"])
assert math.isclose(reported.ProductiveArea, sum(areas.values()))
assert run["ChecksPassed"] and run["SolutionCount"] > 0
hashes_file = output / "source-hashes.json"
if hashes_file.exists():
for name, digest in json.loads(hashes_file.read_text()).items():
assert source_digest(ROOT / name) == digest
print("Saved allocation and layout evidence passed independent checks.")
def write_tables(output):
"""Build the publication tables from the saved evidence, without a solve."""
directory = output / "_tables"
directory.mkdir(exist_ok=True)
def save(name, table, caption):
table = table.copy()
for column in table.columns:
if pd.api.types.is_integer_dtype(table[column]):
table[column] = table[column].map(lambda v: f"{v:,d}")
elif pd.api.types.is_float_dtype(table[column]):
precision = 4 if name == "handling" else 6 if (
column == "Dwell time (years)"
) else 2
table[column] = table[column].map(
lambda v: f"{v:,.{precision}f}"
)
text = table.to_markdown(index=False, disable_numparse=True)
text = text.replace("$", "\\$")
(directory / f"{name}.md").write_text(
text + f"\n\n: {caption} {{#tbl-tut-{name} .striped}}\n",
encoding="utf-8",
)
products, handling, storage, bounds = load_allocation_inputs()
save("products", products.reset_index()[[
"Product", "Demand", "OrderCost", "Price", "SpacePerLoad",
"ReserveShare"
]].rename(columns={
"Demand": "Loads/year", "OrderCost": "Order cost ($)",
"Price": "Price ($/load)", "SpacePerLoad": "m²/load",
"ReserveShare": "Reserve share in RF",
}), "Six-product example inputs. Carrying rate is 10% per year.")
for name, table in (("handling", handling), ("storage", storage)):
save(name, table.rename(columns={str(k): v
for k, v in FLOW_NAMES.items()})
.reset_index(), name.capitalize() + " costs " + (
"($/load/transaction)." if name == "handling"
else "($/average load/year)."
))
save("bounds", bounds.reset_index().rename(columns={
"Lower": "Minimum (m²)", "Upper": "Maximum (m²)"
}), "Functional-area bounds. Total allocated floor area is 100,000 m².")
save("inventory", products.reset_index()[[
"Product", "EOQ", "AverageInventory", "DwellYears"
]].rename(columns={"EOQ": "EOQ (loads)",
"AverageInventory": "Average inventory (loads)",
"DwellYears": "Dwell time (years)"}),
"Calculated order quantities, average inventory, and dwell times.")
assignments = pd.read_csv(output / "assignments.csv")
cost_table = assignments.drop(columns="FlowID")
cost_table.loc[len(cost_table)] = {
"Product": "Total", "Flow": "All",
**{column: assignments[column].sum() for column in
("HandlingCost", "StorageCost", "TotalCost")},
}
save("assignments", cost_table.rename(columns={
"HandlingCost": "Handling ($/year)",
"StorageCost": "Storage ($/year)", "TotalCost": "Total ($/year)"
}), "Optimal product flows and annual cost components for the "
"two-transaction convention.")
save("areas", pd.read_csv(output / "areas.csv").rename(columns={
"Occupied": "Occupied (m²)", "Allocated": "Allocated (m²)",
"UtilizationPct": "Utilization (%)"
}), "One cost-optimal area allocation and its average-inventory use.")
if (output / "departments.csv").exists():
save("departments", pd.read_csv(output / "departments.csv")
.rename(columns={"Area_m2": "Target area (m²)"}),
"Six departments in the reduced Coolblue case.")
runs = pd.read_csv(output / "layout-comparison.csv")
save("layouts", runs[[
"Shape", "Objective", "ProductiveArea", "GrossEnvelope",
"CornerVoid", "Gap"
]].assign(Gap=lambda t: t.Gap.map(lambda g: f"{100*g:.6f}%"))
.rename(columns={
"Objective": "Closeness score (weighted m)",
"ProductiveArea": "Department area (m²)",
"GrossEnvelope": "Envelope (m²)", "CornerVoid": "Void (m²)",
"Gap": "Relative gap"
}), "Computed layouts with the same department areas and weights. "
"The L-case has a different envelope and excluded corner.")
save("diagnostics", runs[[
"Shape", "Status", "Runtime", "Bound", "Gap",
"MaxAreaErrorPct"
]].assign(Gap=lambda t: t.Gap.map(lambda g: f"{100*g:.6f}%"))
.rename(columns={"Runtime": "Runtime (s)",
"Gap": "Relative gap",
"MaxAreaErrorPct": "Max. area error (%)"}),
"Solver diagnostics. Gurobi status 2 means optimal within the "
"requested relative tolerance of 0.0001.")
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--part", choices=["all", "allocation", "layouts"],
default="all")
parser.add_argument("--output", type=Path,
default=ROOT / "tutorial-results")
parser.add_argument("--time-limit", type=int, default=60)
parser.add_argument("--seed", type=int, default=2026)
parser.add_argument("--rf-operations", type=int, choices=[2, 3], default=2)
parser.add_argument("--verify", action="store_true")
args = parser.parse_args()
if args.verify:
verify_saved(args.output, require_layouts=args.part != "allocation")
return
args.output.mkdir(parents=True, exist_ok=True)
if args.part in {"all", "allocation"}:
solve_allocation(args.output, args.rf_operations)
if args.part in {"all", "layouts"}:
solve_layouts(args.output, args.time_limit, args.seed)
if (args.output / "assignments.csv").exists():
write_tables(args.output)
inputs = list((ROOT / "tutorial-inputs").glob("*.csv"))
inputs += list((ROOT / "instances/data/coolblue").glob("*.csv"))
inputs.append(Path(__file__).resolve())
(args.output / "source-hashes.json").write_text(json.dumps({
path.relative_to(ROOT).as_posix():
source_digest(path)
for path in inputs
}, indent=2), encoding="utf-8")
if __name__ == "__main__":
main()