Simultaneous Optimization of Facility Location and Transportation Planning in Supply Chain Management

March 17, 2026

Supply chain design is one of the richest applications of mathematical optimization. A classic — and notoriously challenging — problem is deciding where to open facilities (warehouses, distribution centers) and how to route goods from those facilities to customers, all at the same time. This is known as the Capacitated Facility Location Problem (CFLP) with transportation planning, and it belongs to the NP-hard family of combinatorial optimization problems.

Problem Formulation

Let’s define the problem precisely.

Sets:

$I$ = set of candidate facility sites ${1, \dots, m}$
$J$ = set of customer locations ${1, \dots, n}$

Parameters:

$f_i$ = fixed opening cost of facility $i$
$c_{ij}$ = unit transportation cost from facility $i$ to customer $j$
$d_j$ = demand of customer $j$
$s_i$ = capacity of facility $i$

Decision Variables:

$y_i \in {0, 1}$ — 1 if facility $i$ is opened, 0 otherwise
$x_{ij} \geq 0$ — amount shipped from facility $i$ to customer $j$

Objective (minimize total cost):

$$\min \sum_{i \in I} f_i y_i + \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij}$$

Subject to:

Demand fulfillment for each customer:
$$\sum_{i \in I} x_{ij} = d_j \quad \forall j \in J$$

Capacity constraint per facility:
$$\sum_{j \in J} x_{ij} \leq s_i y_i \quad \forall i \in I$$

Logical constraint (ship only from open facilities):
$$x_{ij} \geq 0, \quad y_i \in {0, 1}$$

This is a Mixed Integer Linear Program (MILP). We solve it with PuLP (CBC solver), then visualize the result in 2D and 3D.

Concrete Example

We set up a realistic scenario:

5 candidate facility sites scattered across a region
12 customers with varying demand levels
Each facility has a different fixed opening cost and capacity
Transportation cost is proportional to Euclidean distance × a cost-per-unit factor

Source Code

# =============================================================
#  Capacitated Facility Location Problem (CFLP)
#  Simultaneous facility opening + transportation planning
#  Solver: PuLP (CBC)  |  Visualization: matplotlib, mpl_toolkits
# =============================================================

# ── 0. Install dependencies (run once in Colab) ──────────────
# !pip install pulp --quiet

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import pulp
import warnings
warnings.filterwarnings("ignore")

np.random.seed(42)

# ══════════════════════════════════════════════════════════════
#  1.  PROBLEM DATA
# ══════════════════════════════════════════════════════════════

# --- Facility candidates (x, y coordinates) ---
facility_coords = np.array([
    [1.5, 8.0],   # F0
    [4.0, 2.5],   # F1
    [7.5, 7.5],   # F2
    [9.0, 3.0],   # F3
    [5.0, 5.5],   # F4
])
facility_names  = ["F0", "F1", "F2", "F3", "F4"]
fixed_costs     = np.array([500, 350, 450, 400, 300])   # opening cost ($)
capacities      = np.array([250, 200, 300, 220, 280])   # max units

# --- Customer locations ---
customer_coords = np.array([
    [0.5, 9.5],  [2.0, 6.5],  [3.5, 9.0],  [1.0, 4.0],
    [4.5, 7.0],  [6.5, 9.0],  [8.5, 8.5],  [7.0, 5.5],
    [9.5, 6.5],  [8.0, 1.5],  [5.5, 3.0],  [3.0, 1.5],
])
customer_names  = [f"C{j}" for j in range(12)]
demands         = np.array([30, 25, 40, 20, 35, 45, 50, 30, 40, 35, 25, 30])

m = len(facility_coords)   # number of facility candidates
n = len(customer_coords)   # number of customers

# --- Transport cost: Euclidean distance × cost rate ---
cost_rate = 10.0   # $ per unit per distance unit
dist = np.linalg.norm(
    facility_coords[:, np.newaxis, :] - customer_coords[np.newaxis, :, :],
    axis=2
)                          # shape (m, n)
transport_cost = cost_rate * dist   # c_ij matrix


# ══════════════════════════════════════════════════════════════
#  2.  MILP MODEL  (PuLP)
# ══════════════════════════════════════════════════════════════

prob = pulp.LpProblem("CFLP", pulp.LpMinimize)

# Decision variables
y = [pulp.LpVariable(f"y_{i}", cat="Binary")        for i in range(m)]
x = [[pulp.LpVariable(f"x_{i}_{j}", lowBound=0)
      for j in range(n)] for i in range(m)]

# Objective: fixed opening costs + transportation costs
prob += (
    pulp.lpSum(fixed_costs[i] * y[i] for i in range(m)) +
    pulp.lpSum(transport_cost[i][j] * x[i][j]
               for i in range(m) for j in range(n))
)

# Constraint 1: each customer's demand must be fully met
for j in range(n):
    prob += pulp.lpSum(x[i][j] for i in range(m)) == demands[j]

# Constraint 2: facility capacity (only if open)
for i in range(m):
    prob += pulp.lpSum(x[i][j] for j in range(n)) <= capacities[i] * y[i]

# Solve (CBC is bundled with PuLP)
solver = pulp.PULP_CBC_CMD(msg=False)
prob.solve(solver)

print(f"Status        : {pulp.LpStatus[prob.status]}")
print(f"Total cost    : ${pulp.value(prob.objective):,.1f}")


# ══════════════════════════════════════════════════════════════
#  3.  EXTRACT RESULTS
# ══════════════════════════════════════════════════════════════

open_flags = [int(round(pulp.value(y[i]))) for i in range(m)]
flow       = np.array([[pulp.value(x[i][j]) for j in range(n)] for i in range(m)])
flow       = np.where(flow is None, 0, flow).astype(float)

open_idx  = [i for i in range(m) if open_flags[i] == 1]
closed_idx = [i for i in range(m) if open_flags[i] == 0]

# Assignment: each customer goes to facility with highest flow
assignment = np.argmax(flow, axis=0)   # shape (n,)

# Per-facility load
loads = [flow[i].sum() for i in range(m)]

# Cost breakdown
fixed_total    = sum(fixed_costs[i] * open_flags[i] for i in range(m))
transport_total = sum(transport_cost[i][j] * flow[i][j]
                      for i in range(m) for j in range(n))

print("\n── Facility decisions ──────────────────────────────────")
for i in range(m):
    status = "OPEN  " if open_flags[i] else "closed"
    print(f"  {facility_names[i]}: {status}  load={loads[i]:.0f}/{capacities[i]}")
print(f"\n  Fixed cost      : ${fixed_total:,.0f}")
print(f"  Transport cost  : ${transport_total:,.0f}")
print(f"  Total cost      : ${fixed_total + transport_total:,.0f}")

# Colour palette (one colour per open facility)
FAC_COLORS = ["#3266ad", "#e05c2e", "#2e9e5b", "#8e44ad", "#c0932f"]

# ══════════════════════════════════════════════════════════════
#  4.  VISUALIZATIONS
# ══════════════════════════════════════════════════════════════

# ---------- helper: arc drawing ----------------------------------
def draw_arc(ax, p0, p1, flow_val, color, lw_max=4):
    """Draw a curved arc proportional to flow between two points."""
    lw = 0.5 + (flow_val / demands.max()) * lw_max
    mid = (p0 + p1) / 2
    perp = np.array([-(p1 - p0)[1], (p1 - p0)[0]]) * 0.08
    ctrl = mid + perp
    # Bezier via many steps
    t = np.linspace(0, 1, 60)
    bx = (1-t)**2*p0[0] + 2*(1-t)*t*ctrl[0] + t**2*p1[0]
    by = (1-t)**2*p0[1] + 2*(1-t)*t*ctrl[1] + t**2*p1[1]
    ax.plot(bx, by, color=color, lw=lw, alpha=0.55, zorder=2)


# ── FIGURE 1: 2-D supply chain network map ──────────────────────
fig1, ax1 = plt.subplots(figsize=(10, 8))
ax1.set_facecolor("#f7f7f5")
fig1.patch.set_facecolor("#f7f7f5")

# Draw flow arcs
for i in open_idx:
    for j in range(n):
        if flow[i][j] > 0.5:
            draw_arc(ax1, facility_coords[i], customer_coords[j],
                     flow[i][j], FAC_COLORS[i])

# Customers
for j in range(n):
    fc = FAC_COLORS[assignment[j]] if assignment[j] in open_idx else "gray"
    ax1.scatter(*customer_coords[j], s=180 + demands[j]*3,
                color=fc, edgecolors="white", linewidths=1.5, zorder=5, marker="o")
    ax1.annotate(customer_names[j], customer_coords[j],
                 textcoords="offset points", xytext=(6, 4),
                 fontsize=8, color="#333")

# Facilities
for i in range(m):
    if open_flags[i]:
        ax1.scatter(*facility_coords[i], s=400, marker="*",
                    color=FAC_COLORS[i], edgecolors="black",
                    linewidths=1.2, zorder=6)
        ax1.annotate(f"{facility_names[i]}\n(load={loads[i]:.0f})",
                     facility_coords[i],
                     textcoords="offset points", xytext=(8, -14),
                     fontsize=9, fontweight="bold", color=FAC_COLORS[i])
    else:
        ax1.scatter(*facility_coords[i], s=200, marker="X",
                    color="#aaaaaa", edgecolors="black",
                    linewidths=0.8, zorder=6, alpha=0.5)
        ax1.annotate(f"{facility_names[i]} (closed)",
                     facility_coords[i],
                     textcoords="offset points", xytext=(6, 4),
                     fontsize=8, color="#aaaaaa")

ax1.set_xlim(-0.5, 10.5)
ax1.set_ylim(-0.5, 11.0)
ax1.set_title(
    f"CFLP Solution: Supply Chain Network Map\n"
    f"Total cost = ${fixed_total+transport_total:,.0f}  "
    f"(fixed ${fixed_total:,} + transport ${transport_total:,.0f})",
    fontsize=12, pad=12
)
ax1.set_xlabel("X coordinate", fontsize=10)
ax1.set_ylabel("Y coordinate", fontsize=10)
ax1.grid(True, color="white", linewidth=0.8)

legend_els = (
    [mlines.Line2D([0],[0], marker="*", color="w", markerfacecolor=FAC_COLORS[i],
                   markersize=12, label=f"{facility_names[i]} (open, cap={capacities[i]})")
     for i in open_idx] +
    [mlines.Line2D([0],[0], marker="X", color="w", markerfacecolor="#aaaaaa",
                   markersize=10, label=f"{facility_names[i]} (closed)")
     for i in closed_idx]
)
ax1.legend(handles=legend_els, loc="lower left", fontsize=8,
           framealpha=0.9, title="Facilities")
plt.tight_layout()
plt.savefig("cflp_network_map.png", dpi=150, bbox_inches="tight")
plt.show()
print("▶ Figure 1 displayed")


# ── FIGURE 2: Cost breakdown bar chart ──────────────────────────
fig2, axes2 = plt.subplots(1, 2, figsize=(12, 5))
fig2.patch.set_facecolor("#f7f7f5")

# Left: stacked bar — fixed vs transport cost per open facility
fac_fix   = [fixed_costs[i] * open_flags[i] for i in range(m)]
fac_trans = [sum(transport_cost[i][j] * flow[i][j] for j in range(n))
             for i in range(m)]
labels_bar = [facility_names[i] for i in range(m)]
x_pos = np.arange(m)

ax = axes2[0]
ax.set_facecolor("#f7f7f5")
bars_f = ax.bar(x_pos, fac_fix,   color=[FAC_COLORS[i] if open_flags[i] else "#cccccc"
                                          for i in range(m)],
                label="Fixed cost", alpha=0.85)
bars_t = ax.bar(x_pos, fac_trans, bottom=fac_fix,
                color=[FAC_COLORS[i] if open_flags[i] else "#eeeeee"
                       for i in range(m)],
                label="Transport cost", alpha=0.5, hatch="//")
ax.set_xticks(x_pos)
ax.set_xticklabels(labels_bar)
ax.set_title("Cost breakdown per facility", fontsize=11)
ax.set_ylabel("Cost ($)")
ax.legend(fontsize=9)
ax.grid(axis="y", color="white", linewidth=0.8)
for i, (f, t) in enumerate(zip(fac_fix, fac_trans)):
    if f + t > 0:
        ax.text(i, f + t + 10, f"${f+t:,.0f}", ha="center", va="bottom", fontsize=8)

# Right: capacity utilization
ax2 = axes2[1]
ax2.set_facecolor("#f7f7f5")
util = [loads[i] / capacities[i] * 100 for i in range(m)]
colors_util = [FAC_COLORS[i] if open_flags[i] else "#cccccc" for i in range(m)]
bars_u = ax2.barh(labels_bar, util, color=colors_util, alpha=0.85, edgecolor="white")
ax2.axvline(100, color="red", lw=1.2, linestyle="--", label="Capacity limit")
ax2.set_xlim(0, 120)
ax2.set_xlabel("Utilization (%)")
ax2.set_title("Facility capacity utilization", fontsize=11)
ax2.legend(fontsize=9)
ax2.grid(axis="x", color="white", linewidth=0.8)
for bar, v in zip(bars_u, util):
    ax2.text(v + 1, bar.get_y() + bar.get_height()/2,
             f"{v:.1f}%", va="center", fontsize=9)

plt.suptitle("CFLP — Cost & Utilization Analysis", fontsize=13, y=1.01)
plt.tight_layout()
plt.savefig("cflp_cost_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
print("▶ Figure 2 displayed")


# ── FIGURE 3: 3-D cost landscape (sensitivity analysis) ─────────
# Vary the cost_rate from 2 to 20 and re-solve to build a surface.
# X-axis: cost_rate | Y-axis: number of open facilities
# Z-axis: total cost

rate_range    = np.linspace(2, 20, 15)
n_open_list   = []
total_cost_list = []

for cr in rate_range:
    tc = cr * dist
    p2 = pulp.LpProblem("CFLP_sweep", pulp.LpMinimize)
    y2 = [pulp.LpVariable(f"y_{i}", cat="Binary") for i in range(m)]
    x2 = [[pulp.LpVariable(f"x_{i}_{j}", lowBound=0) for j in range(n)]
          for i in range(m)]
    p2 += (pulp.lpSum(fixed_costs[i] * y2[i] for i in range(m)) +
           pulp.lpSum(tc[i][j] * x2[i][j]
                      for i in range(m) for j in range(n)))
    for j in range(n):
        p2 += pulp.lpSum(x2[i][j] for i in range(m)) == demands[j]
    for i in range(m):
        p2 += pulp.lpSum(x2[i][j] for j in range(n)) <= capacities[i] * y2[i]
    p2.solve(pulp.PULP_CBC_CMD(msg=False))
    yvals = [int(round(pulp.value(y2[i]))) for i in range(m)]
    n_open_list.append(sum(yvals))
    total_cost_list.append(pulp.value(p2.objective))

# Build grid for surface plot
unique_n = sorted(set(n_open_list))
rate_arr  = np.array(rate_range)
cost_arr  = np.array(total_cost_list)
n_arr     = np.array(n_open_list)

fig3 = plt.figure(figsize=(13, 7))
fig3.patch.set_facecolor("#f7f7f5")

# Left subplot: 3-D scatter + surface interpolation
ax3d = fig3.add_subplot(121, projection="3d")
ax3d.set_facecolor("#f7f7f5")

scatter = ax3d.scatter(rate_arr, n_arr, cost_arr,
                       c=cost_arr, cmap="plasma",
                       s=80, depthshade=True, edgecolors="white", linewidths=0.4)
fig3.colorbar(scatter, ax=ax3d, pad=0.1, shrink=0.6, label="Total cost ($)")

# Highlight current solution point
ax3d.scatter([cost_rate], [sum(open_flags)], [fixed_total + transport_total],
             c="red", s=160, marker="*", zorder=10, label="Current solution")

ax3d.set_xlabel("Transport cost rate", fontsize=9, labelpad=6)
ax3d.set_ylabel("# open facilities", fontsize=9, labelpad=6)
ax3d.set_zlabel("Total cost ($)", fontsize=9, labelpad=6)
ax3d.set_title("3-D Sensitivity:\nCost rate vs # facilities vs Total cost",
               fontsize=10)
ax3d.legend(fontsize=8)

# Right subplot: 3-D bar chart — per-facility transport contribution at baseline
ax3d2 = fig3.add_subplot(122, projection="3d")
ax3d2.set_facecolor("#f7f7f5")

xpos = np.arange(m)           # facility index
ypos = np.arange(n)           # customer index
xposM, yposM = np.meshgrid(xpos, ypos, indexing="ij")
zpos  = np.zeros_like(xposM, dtype=float)
dz    = flow                   # flow[i][j] = height of each bar

dx = dy = 0.6   # bar footprint
cmap_3d = plt.cm.get_cmap("YlOrRd")
max_dz  = dz.max() if dz.max() > 0 else 1.0

for i in range(m):
    for j in range(n):
        if dz[i, j] > 0.1:
            color_val = dz[i, j] / max_dz
            ax3d2.bar3d(xpos[i] - dx/2, ypos[j] - dy/2, 0,
                        dx, dy, dz[i, j],
                        color=cmap_3d(color_val), alpha=0.75,
                        shade=True, edgecolor="white", linewidth=0.2)

ax3d2.set_xticks(xpos)
ax3d2.set_xticklabels(facility_names, fontsize=7)
ax3d2.set_yticks(ypos)
ax3d2.set_yticklabels(customer_names, fontsize=6)
ax3d2.set_xlabel("Facility", fontsize=9, labelpad=6)
ax3d2.set_ylabel("Customer", fontsize=9, labelpad=6)
ax3d2.set_zlabel("Flow (units)", fontsize=9, labelpad=6)
ax3d2.set_title("3-D Flow matrix:\nFacility → Customer shipments",
                fontsize=10)

plt.suptitle("CFLP — 3-D Analysis", fontsize=13)
plt.tight_layout()
plt.savefig("cflp_3d_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
print("▶ Figure 3 displayed")


# ── FIGURE 4: Customer demand vs served cost heatmap ────────────
fig4, ax4 = plt.subplots(figsize=(9, 4))
fig4.patch.set_facecolor("#f7f7f5")
ax4.set_facecolor("#f7f7f5")

# transport cost each customer actually pays
cust_cost = np.array([
    sum(transport_cost[i][j] * flow[i][j] for i in range(m))
    for j in range(n)
])
served_by = [facility_names[assignment[j]] if assignment[j] in open_idx
             else "none" for j in range(n)]

colors_cust = [FAC_COLORS[assignment[j]] if assignment[j] in open_idx
               else "#cccccc" for j in range(n)]

bars = ax4.bar(customer_names, cust_cost, color=colors_cust,
               edgecolor="white", alpha=0.85)
# demand as scatter overlay
ax4_twin = ax4.twinx()
ax4_twin.scatter(customer_names, demands, color="black",
                 zorder=5, s=50, label="Demand (units)")
ax4_twin.set_ylabel("Demand (units)", fontsize=10)
ax4_twin.legend(loc="upper right", fontsize=9)

ax4.set_ylabel("Transport cost ($)", fontsize=10)
ax4.set_title("Per-customer transport cost and demand\n(bar color = serving facility)",
              fontsize=11)
ax4.grid(axis="y", color="white", linewidth=0.8)

legend_els2 = [mpatches.Patch(color=FAC_COLORS[i], label=f"Served by {facility_names[i]}")
               for i in open_idx]
ax4.legend(handles=legend_els2, loc="upper left", fontsize=9)
plt.tight_layout()
plt.savefig("cflp_customer_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
print("▶ Figure 4 displayed")

Code Walkthrough

Section 0 — Dependency

Only PuLP needs to be installed. numpy and matplotlib are pre-installed in Colab. The !pip install pulp --quiet line at the top handles this automatically.

Section 1 — Problem Data

We place 5 facilities and 12 customers on a $[0,10]^2$ map. Each facility has a fixed_cost (paid once if opened) and a capacity ceiling. Each customer has a demand that must be fully satisfied.

The transport cost matrix $c_{ij}$ is computed as:

$$c_{ij} = r \cdot |p_i^{fac} - p_j^{cus}|_2$$

where $r = 10$ is the cost rate per unit per distance unit.

Section 2 — MILP Model

PuLP builds the model symbolically. Binary variables $y_i$ encode open/close decisions; continuous variables $x_{ij}$ encode shipment flows. The lpSum calls translate directly to the mathematical formulation above. The solver is CBC (bundled with PuLP — no external solver needed).

Section 3 — Result Extraction

After solving, we read back variable values with pulp.value(). The argmax trick identifies each customer’s primary serving facility cleanly even when demand is split across multiple facilities.

Section 4 — Visualizations

Figure 1 — Network Map: Bezier-curved arcs represent flows (thicker = larger shipment). Open facilities are shown as stars ★, closed facilities as ✗. Customer dot size scales with demand; dot color matches the serving facility.

Figure 2 — Cost & Utilization: Left panel is a stacked bar (fixed cost base + transport cost overlay) per facility. Right panel is a horizontal bar showing utilization as a percentage of capacity — any bar exceeding 100% would indicate a violated constraint (none should).

Figure 3 — 3-D Analysis (two panels):

Left: A sensitivity sweep over transport cost rate $r \in [2, 20]$. The MILP is re-solved 15 times. The 3-D scatter reveals how the number of open facilities and total cost respond to changing logistics economics. At low cost rates, fewer (larger) facilities are preferred; at high rates, the optimizer opens more facilities to reduce distance.
Right: A 3-D bar chart of the flow matrix $x_{ij}$. Each bar’s height is the volume shipped from facility $i$ to customer $j$. This immediately shows the concentration of supply chains.

Figure 4 — Customer Analysis: Bars show each customer’s actual transport bill; the overlaid scatter shows their demand. Color ties each customer to its serving facility, making service territories visually obvious.

Execution Results

Status        : Optimal
Total cost    : $9,518.1

── Facility decisions ──────────────────────────────────
  F0: OPEN    load=95/250
  F1: OPEN    load=75/200
  F2: OPEN    load=135/300
  F3: OPEN    load=35/220
  F4: OPEN    load=65/280

  Fixed cost      : $2,000
  Transport cost  : $7,518
  Total cost      : $9,518

▶ Figure 1 displayed

▶ Figure 2 displayed

▶ Figure 3 displayed

▶ Figure 4 displayed

Key Insights

The sensitivity analysis in Figure 3 is particularly revealing. At a low transport cost rate, the optimizer is happy to consolidate supply through fewer facilities — paying high transport costs is cheap anyway, so it avoids multiple fixed opening fees. As the transport rate rises, the penalty for long hauls grows, and the solver shifts strategy: open more facilities, place them closer to clusters of demand, and trade fixed costs for transport savings. This crossover behavior is the heart of the simultaneous optimization — you cannot find the right network design without solving both decisions together.