1 Welcome and learning goals

A supply chain is a network of flows: materials flow from suppliers to plants, products flow from plants to distribution centers (DCs), and finished goods flow from DCs to customers. In this lab you will analyze one simulated supply chain with four different tools, because each tool is best at a different job.

Tool	What we use it for in this lab
R	Data simulation, demand forecasting, spatial buffers, network metrics, stress tests
Python	Graph algorithms (`networkx`) and a cost-minimizing transportation model (`scipy`)
Tableau	Interactive maps, flow maps, buffer maps, and a dashboard for decision makers
Sourcemap	Visual, multi-tier supply chain mapping and transparency

By the end of this lab you will be able to:

Compare demand and supply over time and quantify shortfalls and service levels.
Build buffer (service-area) and population-density analyses around facilities.
Describe a supply chain as a network and test how fragile it is.
Extend the analysis to inventory, network design, and emissions.
Move the same data between R, Python, Tableau, and Sourcemap.

Roadmap. We start with the two visual tools (Tableau and Sourcemap) so you can see the finished picture, then build the analysis in R, extend it in Python, and finish with a capstone.

Data note. Everything here is simulated or rounded for teaching. City populations are approximate metro-area figures rounded to one decimal (millions), and the company, plants, and demand are fictional. Do not use these numbers for real decisions.

1.1 Key vocabulary

Term	Plain-English meaning	Example from our case
Node	A place where something is made, stored, or consumed	A plant in Columbus
Edge	A flow between two nodes	Columbus plant to Atlanta DC
Tier	A layer of the chain	Supplier, Plant, DC, Customer
Lead time	Time between ordering and receiving	2 weeks from plant to DC
Fill rate	Share of demand actually delivered	11,000 delivered of 12,000 ordered = 91.7%
Buffer	The area within a set distance of a facility	Everything within 300 miles of Denver
Centrality	How important a node is in the network	A plant that most paths pass through

Code is folded by default: click the Code button above any output to reveal the R or Python syntax.

1.2 Packages

Run this once (it is not evaluated when knitting).

install.packages(c(
  "tidyverse", "sf", "maps", "tidygraph", "ggraph", "igraph",
  "kableExtra", "scales", "lubridate", "checkdown", "reticulate"
))

# Python side (only needed for the Python section)
reticulate::py_install(c("pandas", "numpy", "networkx", "scipy"))

library(tidyverse)
library(sf)
library(tidygraph)
library(ggraph)
library(igraph)
library(kableExtra)
library(scales)
library(lubridate)
library(checkdown)
library(reticulate)

# Quiz helper (wraps checkdown). checkdown pastes the answer and feedback text into
# JavaScript strings wrapped in single quotes, so a straight apostrophe breaks the
# quiz (the Check button then loads a blank page). We swap them for curly ones.
quiz <- function(answer, options, right = "Correct!", wrong = "Not quite. Try again.") {
  clean <- function(x) gsub("'", "\u2019", x, fixed = TRUE)
  check_question(
    answer = clean(answer), options = clean(options), type = "radio",
    alignment = "vertical", right = clean(right), wrong = clean(wrong)
  )
}

# A small helper so every table in the lab looks the same
nice_table <- function(df, caption = NULL, digits = 1) {
  df %>%
    kbl(caption = caption, digits = digits, format.args = list(big.mark = ",")) %>%
    kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
}

2 Tableau: interactive maps and a decision dashboard

Where does the data come from? The CSV files used in the next two sections (nodes.csv, edges.csv, national_weekly.csv, and so on) are created by the export chunk at the end of the R analysis (“Export data for Python, Tableau, and Sourcemap”). Either work through the R sections first and come back here, or use the data/ folder provided by your instructor. You can start exploring the tools now, because the steps do not depend on the numbers.

Tableau cannot be run from R Markdown, so this section is a guided build. Use Tableau Public (free) or Tableau Desktop and connect to the CSV files in data/ (Connect → Text file).

2.1 Build 1: Demand vs. supply (use `national_weekly.csv`)

Concept: dual axis. Two measures drawn on one shared axis so you can compare them directly. Example: demand (blue) and supply (red) on the same weekly timeline. Whenever blue rises above red, you have a shortfall.

Concept: calculated field. A new column defined by a formula. Example: the Status field below labels each week “Shortfall” or “OK” so you can color the chart by it.

Drag week to Columns (continuous, by Week).
Drag demand and supply to Rows, then make a dual axis (right-click one axis → Dual Axis → Synchronize Axis).
Create a calculated field Status:

IF [supply] < [demand] THEN "Shortfall" ELSE "OK" END

Put Status on Color. Add a reference band over the disruption weeks.

2.2 Build 2: Network map (use `nodes.csv`)

Concept: geographic role. A label that tells Tableau a number is a latitude or longitude so it can place a point on a map. Example: once lat and lon have geographic roles, each row of nodes.csv becomes a dot on the map.

Double-click lon and lat (set geographic roles to Longitude and Latitude).
Drag type to Color and pop_m to Size; use Shape to make DCs triangles.

2.3 Build 3: Flow map (use `edges.csv`)

Concept: flow map. A map where each row of data is a link between two places, drawn as a line. Example: the row “Columbus plant to Atlanta DC” becomes a line from Columbus to Atlanta. Thick, busy corridors show where the chain depends on a few routes.

Create a line from each origin to its destination:

MAKELINE(MAKEPOINT([from_lat], [from_lon]), MAKEPOINT([to_lat], [to_lon]))

Drag the new field to Detail, set the mark type to Line, and color by tier. Layer it with Build 2 as a dual-axis map.

2.4 Build 4: Interactive buffer map

Concept: parameter. A value the viewer can change with a slider or box. Example: slide Radius (mi) from 300 to 500 and every DC circle grows instantly, so you can see which customers fall inside.

Create a parameter Radius (mi) (integer, 50–600, step 50) and show the control.
Create a calculated field (filtered to type = "DC"):

BUFFER(MAKEPOINT([lat], [lon]), [Radius (mi)], 'mi')

Drag it to Detail on a map layer, lower the opacity, and layer the customers on top (Build 2).

2.5 Build 5: Dashboard

Concept: filter action. Clicking a mark on one sheet filters the other sheets. Example: click the Chicago DC on the map and the demand chart shows only the customers that Chicago serves.

Combine Builds 1–4 on one dashboard. Add a Filter action (click a DC, filter the demand view) and a tooltip with pop_m and distance.

Tableau + R/Python. Tableau Desktop can call R or Python through external services (Rserve or TabPy) for table calculations. This is not available in Tableau Public, which is why this lab passes data through CSV files.

Quick check. Which Tableau function draws a circle of a given radius around a point?

3 Sourcemap: multi-tier supply chain mapping

Our R network has three tiers. Real supply chains go much deeper (components, raw materials, farms, mines), and that depth is where Sourcemap is useful: it is built for mapping and communicating multi-tier supply chains visually.

Concept: tiers and traceability. Tier 1 suppliers sell directly to you; tier 2 suppliers sell to your tier 1 suppliers, and so on back to raw materials. Traceability means being able to follow a product back through those tiers. Example (chocolate): the chocolate maker’s tier 1 supplier is a cocoa processor, tier 2 is the cooperative that collects beans, and tier 3 is the individual farms. Most companies can name tier 1 but know little about tiers 2 and 3, which is where hidden risk lives.

Sign in to Sourcemap and create a new supply chain map.
Open the import option and look at the current import template (templates change, so use the live one, not a screenshot).
Map the columns in data/sourcemap_nodes.csv (Name, Type, Latitude, Longitude) and data/sourcemap_edges.csv (From, To, Tier, Miles) onto the template columns, and import.
Compare the Sourcemap view with your R network plot and your Tableau flow map.
Add at least one tier-2 supplier for each supplier (e.g., a raw-material source) and note what new risks appear.

Reflection questions

What does Sourcemap show more easily than R or Tableau? What does it show less easily?
Which questions (forecasting, optimization, stress tests) still require R or Python?

4 Our case: a simulated national supply chain

Concept: simulated data. Real company data is confidential, so we generate realistic data with random numbers. set.seed() fixes the random draws so everyone gets identical results. Example: with set.seed(580), every student’s “random” demand series is the same, which makes answers comparable.

Our fictional company buys components from 4 suppliers, assembles them in 2 plants, stores them in 5 DCs, and sells to 20 metro-area customers.

suppliers <- tribble(
  ~id,  ~label,    ~lon,    ~lat,
  "S1", "Houston", -95.37,  29.76,
  "S2", "Detroit", -83.05,  42.33,
  "S3", "Charlotte", -80.84, 35.23,
  "S4", "Phoenix", -112.07, 33.45
) %>% mutate(type = "Supplier")

plants <- tribble(
  ~id,  ~label,    ~lon,   ~lat,  ~capacity,
  "P1", "Columbus", -83.00, 39.96, 9000,
  "P2", "Dallas",   -96.80, 32.78, 9000
) %>% mutate(type = "Plant")

dcs <- tribble(
  ~id,  ~label,    ~lon,    ~lat,
  "D1", "Atlanta", -84.39,  33.75,
  "D2", "Chicago", -87.63,  41.88,
  "D3", "Ontario CA", -117.65, 34.06,
  "D4", "Denver",  -104.99, 39.74,
  "D5", "Newark",  -74.17,  40.74
) %>% mutate(type = "DC")

# pop_m = approximate metro population in millions (rounded, for teaching)
customers <- tribble(
  ~label,           ~lon,     ~lat,   ~pop_m,
  "New York",       -74.00,   40.71,  19.5,
  "Los Angeles",    -118.24,  34.05,  12.8,
  "Chicago",        -87.63,   41.88,   9.4,
  "Dallas",         -96.80,   32.78,   8.1,
  "Houston",        -95.37,   29.76,   7.5,
  "Washington DC",  -77.04,   38.91,   6.3,
  "Philadelphia",   -75.17,   39.95,   6.2,
  "Atlanta",        -84.39,   33.75,   6.2,
  "Miami",          -80.19,   25.76,   6.1,
  "Phoenix",        -112.07,  33.45,   5.0,
  "Boston",         -71.06,   42.36,   4.9,
  "San Francisco",  -122.42,  37.77,   4.6,
  "Detroit",        -83.05,   42.33,   4.4,
  "Seattle",        -122.33,  47.61,   4.0,
  "Minneapolis",    -93.27,   44.98,   3.7,
  "San Diego",      -117.16,  32.72,   3.3,
  "Denver",         -104.99,  39.74,   3.0,
  "St. Louis",      -90.20,   38.63,   2.8,
  "Charlotte",      -80.84,   35.23,   2.8,
  "Kansas City",    -94.58,   39.10,   2.2
) %>%
  mutate(id = sprintf("C%02d", row_number()), type = "Customer")

nodes <- bind_rows(suppliers, plants, dcs, customers) %>%
  select(id, label, type, lon, lat, pop_m, capacity)

nodes %>% count(type) %>% nice_table("Nodes in our network", digits = 0)

Nodes in our network
type	n
Customer	20
DC	5
Plant	2
Supplier	4

4.1 Spatial objects (`sf`)

Concept: coordinate reference system (CRS). A rule for turning positions on a round Earth into numbers. Latitude/longitude (EPSG:4326) is measured in degrees, and a degree is not a fixed distance. Example: one degree of longitude is about 69 miles at the equator but only about 53 miles at 40 degrees north. A projected CRS such as EPSG:5070 uses meters, so distances and areas are reliable.

Distances and areas must be computed in a projected coordinate system (meters), not in raw latitude/longitude degrees. We use EPSG:5070 (NAD83 / Conus Albers), an equal-area projection for the contiguous U.S.

crs_proj <- 5070
mi <- 1609.344                       # meters per mile

nodes_sf   <- st_as_sf(nodes, coords = c("lon", "lat"), crs = 4326, remove = FALSE)
nodes_proj <- st_transform(nodes_sf, crs_proj)

suppliers_p <- filter(nodes_proj, type == "Supplier")
plants_p    <- filter(nodes_proj, type == "Plant")
dcs_p       <- filter(nodes_proj, type == "DC")
cust_p      <- filter(nodes_proj, type == "Customer")

total_pop <- sum(customers$pop_m)

Quick check. Why do we transform to a projected CRS before calling st_buffer() and st_area()?

5 Demand and supply analysis

5.1 Simulate three years of weekly demand

Concept: trend, seasonality, and noise. Most demand series combine three parts: a trend (long-run growth or decline), seasonality (a repeating yearly pattern), and noise (random variation). Example (toys): trend = a growing customer base, seasonality = the December spike, noise = a random slow week in March.

Demand for each customer is proportional to its population, then shaped by a trend, annual seasonality, a holiday bump (ISO weeks 47–51), and random noise.

set.seed(580)
weeks <- seq(as.Date("2023-01-02"), by = "week", length.out = 156)

demand_weekly <- customers %>%
  select(id, label, pop_m) %>%
  tidyr::crossing(week = weeks) %>%
  mutate(
    t        = as.numeric(week - min(weeks)) / 7,
    woy      = isoweek(week),
    trend    = 1 + 0.0015 * t,
    season   = 1 + 0.25 * sin(2 * pi * t / 52),
    holiday  = if_else(woy %in% 47:51, 0.15, 0),
    noise    = rlnorm(n(), meanlog = 0, sdlog = 0.12),
    demand   = round(pop_m * 90 * trend * (season + holiday) * noise)
  )

demand_weekly %>% slice_head(n = 5) %>% select(id, label, week, demand) %>% nice_table("First rows of demand data", 0)

First rows of demand data
id	label	week	demand
C01	New York	2023-01-02	1,769
C01	New York	2023-01-09	2,025
C01	New York	2023-01-16	1,732
C01	New York	2023-01-23	1,686
C01	New York	2023-01-30	1,589

5.2 Simulate supply and a disruption

Concept: capacity and disruption. Capacity is the most a facility can produce per period. A disruption temporarily cuts it. Example: a plant that normally makes 9,000 units a week loses 70% of its output in a fire, so it makes only 2,700 for six weeks.

Each plant can make 9,000 units a week. In weeks 80–85 Plant P2 (Dallas) loses 70% of its output (think: a fire or a labor strike).

supply_weekly <- tibble(week = weeks, t = seq_along(weeks)) %>%
  mutate(
    supply = sum(plants$capacity) -
      if_else(between(t, 80, 85), 0.70 * 9000, 0)
  ) %>%
  select(week, supply)

national <- demand_weekly %>%
  group_by(week) %>%
  summarise(demand = sum(demand), .groups = "drop") %>%
  left_join(supply_weekly, by = "week") %>%
  mutate(net = supply - demand)

5.3 Inventory, shortfalls, and service level

Concept: inventory balance and fill rate. Each week: ending inventory = starting inventory + supply - demand (never below zero). Demand you cannot cover is unmet (lost sales). Fill rate = 1 - unmet / demand. Example: demand 12,000, supply 10,000, starting inventory 1,000. Inventory covers 1,000 of the 2,000 gap, so unmet = 1,000 and fill rate = 1 - 1,000 / 12,000 = 91.7%.

Surplus production builds inventory; shortfalls are first covered from inventory and any remainder is lost (unmet demand). Because each week’s inventory depends on the previous week, we use purrr::accumulate() rather than a loop.

inv_end <- accumulate(national$net, ~ pmax(0, .x + .y), .init = 0)[-1]

national <- national %>%
  mutate(
    inv_end   = inv_end,
    inv_start = lag(inv_end, default = 0),
    unmet     = pmax(0, demand - supply - inv_start),
    service   = 1 - unmet / demand
  )

tibble(
  Metric = c("Average weekly demand", "Average weekly supply",
             "Weeks with unmet demand", "Overall fill rate"),
  Value  = c(mean(national$demand), mean(national$supply),
             sum(national$unmet > 0), 1 - sum(national$unmet) / sum(national$demand))
) %>%
  mutate(Value = c(comma(Value[1:3], accuracy = 1), percent(Value[4], accuracy = 0.1))) %>%
  nice_table("Supply chain KPIs")

Supply chain KPIs
Metric	Value
Average weekly demand	12,568
Average weekly supply	17,758
Weeks with unmet demand	0
Overall fill rate	100.0%

national %>%
  select(week, Demand = demand, Supply = supply) %>%
  pivot_longer(-week, names_to = "series", values_to = "units") %>%
  ggplot(aes(week, units, color = series)) +
  annotate("rect", xmin = weeks[80], xmax = weeks[85], ymin = -Inf, ymax = Inf,
           alpha = 0.15, fill = "red") +
  geom_line(linewidth = 0.9) +
  scale_y_continuous(labels = comma) +
  scale_color_manual(values = c(Demand = "#1f77b4", Supply = "#d62728")) +
  labs(title = "Weekly demand vs. supply capacity",
       subtitle = "Shaded band = Plant P2 disruption",
       x = NULL, y = "Units per week", color = NULL) +
  theme_minimal(base_size = 13)

Quick check. Supply capacity is above average demand, yet we still lose sales. Which explanation fits best?

5.4 Forecast the next quarter

Concept: Holt-Winters forecasting. A method that keeps running estimates of the level, trend, and seasonal pattern, and updates them as each week arrives. A prediction interval is the range where future demand should fall most of the time. Example: a forecast of 14,000 with a 95% interval of 12,000 to 16,000 means we expect demand in that range about 95 weeks out of 100.

We use Holt-Winters exponential smoothing from base R (no extra packages) with multiplicative seasonality.

demand_ts <- ts(national$demand, frequency = 52)
hw        <- HoltWinters(demand_ts, seasonal = "multiplicative")
fc        <- predict(hw, n.ahead = 13, prediction.interval = TRUE)

fc_df <- tibble(
  week     = seq(max(weeks) + 7, by = "week", length.out = 13),
  forecast = as.numeric(fc[, "fit"]),
  lwr      = as.numeric(fc[, "lwr"]),
  upr      = as.numeric(fc[, "upr"])
)

capacity <- sum(plants$capacity)

ggplot() +
  geom_line(data = tail(national, 52), aes(week, demand), color = "grey30") +
  geom_ribbon(data = fc_df, aes(week, ymin = lwr, ymax = upr), fill = "#1f77b4", alpha = 0.2) +
  geom_line(data = fc_df, aes(week, forecast), color = "#1f77b4", linewidth = 1) +
  geom_hline(yintercept = capacity, linetype = "dashed", color = "#d62728") +
  scale_y_continuous(labels = comma) +
  labs(title = "13-week demand forecast vs. plant capacity",
       subtitle = "Dashed line = total plant capacity",
       x = NULL, y = "Units per week") +
  theme_minimal(base_size = 13)

fc_df %>% mutate(over_capacity_risk = upr > capacity) %>% count(over_capacity_risk)

## # A tibble: 1 × 2
##   over_capacity_risk     n
##   <lgl>              <int>
## 1 FALSE                 13

Discussion. In how many forecast weeks does the upper prediction bound exceed capacity? What would you recommend: build inventory, add capacity, or accept lost sales?

6 Buffer and population-density analysis

A buffer is the area within a given distance of a facility. Buffers answer questions like “How many people can a DC reach within a one-day truck drive?”

6.1 Buffers around DCs

Concept: buffer (service area). All locations within a chosen distance of a facility. Example: 300 miles is roughly a one-day truck trip, so a 300-mile buffer around the Denver DC shows who can get next-day delivery.

Concept: diminishing returns. Area grows with the square of the radius (\(\pi r^2\)), but people are concentrated in cities, so extra distance mostly adds empty land. Example: doubling the radius from 150 to 300 miles grows the area from about 70,700 to about 282,700 square miles (four times larger), but the number of new customers is far smaller than four times.

radii <- c(150, 300, 500)   # miles

# Union of all DC buffers = total network coverage
coverage <- map_dfr(radii, function(r) {
  buf     <- st_union(st_buffer(dcs_p, dist = r * mi))
  covered <- lengths(st_intersects(cust_p, buf)) > 0
  tibble(
    radius_mi         = r,
    customers_covered = sum(covered),
    pop_covered_m     = sum(cust_p$pop_m[covered]),
    pct_pop_covered   = pop_covered_m / total_pop,
    buffer_area_sqmi  = as.numeric(st_area(buf)) / mi^2
  )
}) %>%
  mutate(covered_pop_per_sqmi = pop_covered_m * 1e6 / buffer_area_sqmi)

coverage %>%
  mutate(pct_pop_covered = percent(pct_pop_covered, accuracy = 0.1)) %>%
  nice_table("Network coverage by buffer radius", digits = 1)

Network coverage by buffer radius
radius_mi	customers_covered	pop_covered_m	pct_pop_covered	buffer_area_sqmi	covered_pop_per_sqmi
150	7	60.4	49.2%	353,267.7	171.0
300	12	81.6	66.4%	1,412,587.3	57.8
500	16	97.1	79.1%	3,370,864.7	28.8

Notice the diminishing returns: each larger radius covers more area, but the number of additional people shrinks. covered_pop_per_sqmi is a coverage-efficiency measure (customers served per square mile of buffer), not a census density.

6.2 Reach of each DC

Concept: reach. The population one DC can serve inside a given radius, on its own. Comparing DCs shows which locations earn their keep. Example: if Newark reaches 40 million people within 300 miles and Denver reaches 6 million, Newark is serving a far denser market.

dc_reach <- map_dfr(radii, function(r) {
  hits <- st_is_within_distance(dcs_p, cust_p, dist = r * mi)
  tibble(DC = dcs_p$label, radius_mi = r,
         pop_reached_m = map_dbl(hits, ~ sum(cust_p$pop_m[.x])))
})

dc_reach %>%
  pivot_wider(names_from = radius_mi, values_from = pop_reached_m, names_prefix = "within_") %>%
  nice_table("Population (millions) reachable from each DC")

Population (millions) reachable from each DC
DC	within_150	within_300	within_500
Atlanta	6.2	9.0	11.8
Chicago	9.4	16.6	22.5
Ontario CA	16.1	16.1	25.7
Denver	3.0	3.0	3.0
Newark	25.7	36.9	41.3

6.3 Map the buffers and find underserved customers

Concept: underserved customers. Customers outside every buffer. Example: a city 450 miles from the nearest DC cannot be reached in one day, so it may need a new DC, a partner warehouse, or a slower delivery promise.

states <- st_as_sf(maps::map("state", fill = TRUE, plot = FALSE))

buf300 <- st_buffer(dcs_p, 300 * mi) %>% st_transform(4326)

# Distance (miles) from each customer to its nearest DC
cust_p <- cust_p %>%
  mutate(nearest_dc_mi = apply(st_distance(., dcs_p), 1, min) / mi,
         served_300 = nearest_dc_mi <= 300)

cust_plot <- st_transform(cust_p, 4326)
dcs_plot  <- st_transform(dcs_p, 4326)

ggplot() +
  geom_sf(data = states, fill = "grey97", color = "grey70") +
  geom_sf(data = buf300, fill = "#1f77b4", alpha = 0.12, color = "#1f77b4") +
  geom_sf(data = cust_plot, aes(size = pop_m, color = served_300), alpha = 0.8) +
  geom_sf(data = dcs_plot, shape = 17, size = 4, color = "black") +
  scale_color_manual(values = c(`TRUE` = "#2ca02c", `FALSE` = "#d62728"),
                     labels = c(`TRUE` = "Within 300 mi", `FALSE` = "Beyond 300 mi")) +
  scale_size_continuous(range = c(2, 9), name = "Pop. (millions)") +
  labs(title = "300-mile DC service areas", color = NULL) +
  theme_void(base_size = 13)

cust_p %>% st_drop_geometry() %>%
  filter(!served_300) %>% select(label, pop_m, nearest_dc_mi) %>%
  arrange(desc(pop_m)) %>% nice_table("Customers beyond 300 miles of any DC")

Customers beyond 300 miles of any DC
label	pop_m	nearest_dc_mi
Dallas	8.1	662.6
Houston	7.5	699.8
Miami	6.1	605.2
Phoenix	5.0	321.7
San Francisco	4.6	369.9
Seattle	4.0	971.1
Minneapolis	3.7	354.9
Kansas City	2.2	411.2

6.4 Optional: real population density with Census data

Concept: population density. People per square mile (density = population / area). Example: a county with 100,000 people in 500 square miles has a density of 200 people per square mile. Dense areas need frequent deliveries to many addresses; sparse areas need long, less efficient routes.

For a true density analysis (people per square mile by county), use tidycensus. This needs a free Census API key, so it is not evaluated here.

library(tidycensus)
census_api_key("YOUR_KEY_HERE", install = TRUE)

county_pop <- get_acs(
  geography = "county", variables = "B01003_001",
  year = 2022, survey = "acs5", geometry = TRUE
) %>%
  st_transform(crs_proj) %>%
  mutate(area_sqmi = as.numeric(st_area(.)) / mi^2,
         density = estimate / area_sqmi)

# How many people (and what average density) live inside the 300-mile buffer?
buf_union <- st_union(st_buffer(dcs_p, 300 * mi))
inside <- county_pop[lengths(st_intersects(st_centroid(county_pop), buf_union)) > 0, ]
sum(inside$estimate)

Quick check. Going from a 300-mile to a 500-mile buffer increases the buffer area a lot but adds relatively few customers. What does this tell a network planner?

7 Network analysis

7.1 Build the network

Concept: nodes, edges, and direction. A node is a facility, an edge is a shipping link, and a directed edge has a from-to orientation, which for us follows the flow of goods. Redundancy means a node has more than one source. Example: each customer is linked to a primary DC (rank 1) and a backup DC (rank 2), so one DC failure does not strand it.

We connect the tiers using k-nearest links to give the network some redundancy: each plant is supplied by its 3 nearest suppliers, each DC by its 2 nearest plants, and each customer by its 2 nearest DCs (rank 1 = primary, rank 2 = backup).

link_nearest <- function(src, dst, k) {
  d <- st_distance(dst, src)
  d <- matrix(as.numeric(d), nrow = nrow(d))   # rows = dst, cols = src (meters)
  map_dfr(seq_len(nrow(d)), function(i) {
    idx <- order(d[i, ])[seq_len(k)]
    tibble(from  = src$id[idx],
           to    = dst$id[i],
           # 1-mile floor: a DC and customer in the same city would otherwise have
           # a 0-mile edge, and igraph weights must be strictly positive
           miles = pmax(d[i, idx] / mi, 1),
           rank  = seq_len(k))
  })
}

edges <- bind_rows(
  link_nearest(suppliers_p, plants_p, k = 3) %>% mutate(tier = "Supplier to Plant"),
  link_nearest(plants_p,    dcs_p,    k = 2) %>% mutate(tier = "Plant to DC"),
  link_nearest(dcs_p,       cust_p,   k = 2) %>% mutate(tier = "DC to Customer")
)

g <- tbl_graph(
  nodes = nodes %>% select(id, label, type, lon, lat, pop_m),
  edges = edges, directed = TRUE, node_key = "id"
)

edges %>% count(tier) %>% nice_table("Edges by tier", 0)

Edges by tier
tier	n
DC to Customer	40
Plant to DC	10
Supplier to Plant	6

ggraph(g, layout = "manual", x = lon, y = lat) +
  geom_edge_link(aes(colour = tier), alpha = 0.45) +
  geom_node_point(aes(colour = type), size = 3.5) +
  geom_node_text(aes(label = ifelse(type != "Customer", label, "")),
                 repel = TRUE, size = 3.3) +
  labs(title = "Supply chain network (positions = real coordinates)",
       colour = NULL, edge_colour = NULL) +
  theme_void(base_size = 13)

7.2 Which nodes matter most?

Concept: degree. The number of links a node has. In-degree counts links arriving (sources), out-degree counts links leaving (destinations). Example: a DC with in-degree 2 can be supplied by two plants.

Concept: betweenness. How many shortest supplier-to-customer routes pass through a node. Example: if nearly every route from western suppliers to eastern customers runs through one plant, that plant has very high betweenness and is a likely single point of failure.

Centrality metrics rank nodes by structural importance. Betweenness counts how many shortest paths (weighted by miles) pass through a node, so a high value means “many flows depend on me”.

g <- g %>%
  mutate(
    in_degree  = centrality_degree(mode = "in"),
    out_degree = centrality_degree(mode = "out"),
    betweenness = centrality_betweenness(weights = miles, directed = TRUE)
  )

g %>% as_tibble() %>%
  filter(type %in% c("Supplier", "Plant", "DC")) %>%
  select(label, type, in_degree, out_degree, betweenness) %>%
  arrange(desc(betweenness)) %>%
  nice_table("Centrality of supply chain facilities")

Centrality of supply chain facilities
label	type	in_degree	out_degree	betweenness
Columbus	Plant	3	5	51
Dallas	Plant	3	5	49
Atlanta	DC	2	8	30
Denver	DC	2	10	29
Chicago	DC	2	9	27
Newark	DC	2	7	21
Ontario CA	DC	2	6	13
Houston	Supplier	0	2	0
Detroit	Supplier	0	1	0
Charlotte	Supplier	0	2	0
Phoenix	Supplier	0	1	0

Quick check. A plant with very high betweenness is best described as:

7.3 Stress test: what if facilities fail?

Concept: stress test (reachability). Remove nodes from the network and check who can still be reached from a supplier. Example: if both plants (Columbus and Dallas) fail, no supplier can reach any DC, so 0% of customers are served. A pair like two neighboring DCs may strand only the customers that depended on both.

We remove every single facility, then every pair of facilities, and measure the share of customer population that can still be reached from at least one supplier.

ig <- as.igraph(g)

pct_served <- function(graph) {
  src  <- V(graph)[type == "Supplier"]
  cust <- V(graph)[type == "Customer"]
  d    <- distances(graph, v = src, to = cust, mode = "out")
  ok   <- colSums(is.finite(d)) > 0
  sum(cust$pop_m[ok]) / total_pop
}

labs_lookup <- setNames(V(ig)$label, V(ig)$id)
infra_ids   <- V(ig)$id[V(ig)$type != "Customer"]

stress_test <- function(m) {
  combn(infra_ids, m, simplify = FALSE) %>%
    map_dfr(function(ids) {
      g2 <- delete_vertices(ig, which(V(ig)$id %in% ids))
      tibble(scenario = paste(labs_lookup[ids], collapse = " + "),
             failures = m, pct_served = pct_served(g2))
    })
}

stress <- bind_rows(stress_test(1), stress_test(2))

stress %>%
  group_by(failures) %>% slice_min(pct_served, n = 5, with_ties = FALSE) %>%
  ungroup() %>%
  mutate(pct_served = percent(pct_served, accuracy = 0.1)) %>%
  nice_table("Worst-case failure scenarios", 0)

Worst-case failure scenarios
scenario	failures	pct_served
Houston	1	100.0%
Detroit	1	100.0%
Charlotte	1	100.0%
Phoenix	1	100.0%
Columbus	1	100.0%
Columbus + Dallas	2	0.0%
Chicago + Newark	2	71.5%
Ontario CA + Denver	2	73.4%
Atlanta + Chicago	2	85.0%
Atlanta + Denver	2	87.3%

Discussion. Which single facility is most critical? Which pair is worst? How would you add redundancy (a third DC link per customer? a third plant?) and what would it cost?

8 Beyond the basics: inventory, network design, and emissions

8.1 Safety stock at each DC

Concept: safety stock and reorder point. Safety stock is extra inventory held to survive demand surprises while waiting for a replenishment. The reorder point is the inventory level that triggers a new order: average demand during the lead time plus safety stock. Example: weekly demand averages 5,000 units (standard deviation 1,000), lead time is 4 weeks, and \(z = 1.645\) (95% service). Safety stock = \(1.645 \times 1{,}000 \times \sqrt{4}\) = 3,290 units. Reorder point = \(5{,}000 \times 4 + 3{,}290\) = 23,290 units.

Safety stock protects against demand variability during the replenishment lead time:

\[SS = z \times \sigma_{weekly} \times \sqrt{L}\]

where \(z\) is the service-level factor (95% gives \(z \approx 1.645\)), \(\sigma_{weekly}\) is the standard deviation of weekly demand, and \(L\) is the lead time in weeks. We assume lead time = 1 week + 1 week per 1,000 miles from the primary plant.

primary_dc <- edges %>% filter(tier == "DC to Customer", rank == 1) %>%
  select(dc = from, customer = to)

dc_stats <- demand_weekly %>%
  left_join(primary_dc, by = c("id" = "customer")) %>%
  group_by(dc, week) %>% summarise(demand = sum(demand), .groups = "drop") %>%
  group_by(dc) %>% summarise(mean_wk = mean(demand), sd_wk = sd(demand), .groups = "drop")

dc_lead <- edges %>% filter(tier == "Plant to DC", rank == 1) %>%
  select(dc = to, plant = from, plant_miles = miles)

safety <- dc_stats %>%
  left_join(dc_lead, by = "dc") %>%
  left_join(select(nodes, id, DC = label), by = c("dc" = "id")) %>%
  mutate(lead_wks      = 1 + plant_miles / 1000,
         safety_stock  = qnorm(0.95) * sd_wk * sqrt(lead_wks),
         reorder_point = mean_wk * lead_wks + safety_stock)

safety %>% select(DC, mean_wk, sd_wk, lead_wks, safety_stock, reorder_point) %>%
  nice_table("Safety stock and reorder point by DC (95% service level)")

Safety stock and reorder point by DC (95% service level)
DC	mean_wk	sd_wk	lead_wks	safety_stock	reorder_point
Atlanta	2,307.5	418.1	1.4	824.8	4,144.8
Chicago	2,299.2	411.3	1.3	764.0	3,695.3
Ontario CA	3,034.8	567.3	2.2	1,383.2	8,052.6
Denver	1,132.0	207.4	1.7	440.0	2,322.1
Newark	3,794.0	683.1	1.5	1,359.9	6,918.4

Quick check. Our \(\sigma_{weekly}\) includes seasonal swings, not just random noise. What is the likely effect on safety stock?

8.2 Network design: how many DCs do we really need?

Concept: facility location (p-median). Choose p facilities so the demand-weighted distance to the nearest one is as small as possible. Example: with 5 candidate DCs, we test every possible 1-DC, 2-DC, 3-DC, 4-DC, and 5-DC network and keep the best of each size. More DCs cut shipping distance but add fixed costs.

Concept: ton-miles and emissions. Ton-miles = tons shipped x miles. Emissions = ton-miles x an emission factor. Example: 10 tons shipped 500 miles = 5,000 ton-miles; at a placeholder factor of 0.1 kg CO2e per ton-mile, that is 500 kg CO2e.

We test every combination of DCs (a brute-force p-median problem) and keep the best set for each network size, measured by demand-weighted miles to the nearest DC. We then convert miles into an emissions estimate.

# Parameters you can change (assumptions, not facts!)
tons_per_unit     <- 0.02    # shipping weight of one unit (tons)
ef_kg_per_tonmile <- 0.10    # PLACEHOLDER emission factor (kg CO2e per ton-mile);
                             # replace with a current EPA or GLEC value

cust_avg <- demand_weekly %>% group_by(id) %>% summarise(avg_demand = mean(demand), .groups = "drop")
w <- cust_avg$avg_demand[match(cust_p$id, cust_avg$id)]

d_cd <- st_distance(cust_p, dcs_p)
d_cd <- matrix(as.numeric(d_cd) / mi, nrow = nrow(d_cd), dimnames = list(cust_p$id, dcs_p$id))

pmedian <- map_dfr(1:5, function(k) {
  map_dfr(combn(colnames(d_cd), k, simplify = FALSE), function(set) {
    nearest <- apply(d_cd[, set, drop = FALSE], 1, min)
    tibble(n_dcs = k, dcs = paste(labs_lookup[set], collapse = ", "),
           wtd_miles = sum(w * nearest) / sum(w))
  }) %>% slice_min(wtd_miles, n = 1, with_ties = FALSE)
}) %>%
  mutate(weekly_tonmiles = wtd_miles * sum(w) * tons_per_unit,
         weekly_co2e_t   = weekly_tonmiles * ef_kg_per_tonmile / 1000)

pmedian %>% nice_table("Best DC set for each network size")

Best DC set for each network size
n_dcs	dcs	wtd_miles	weekly_tonmiles	weekly_co2e_t
1	Chicago	888.8	223,410.9	22.3
2	Ontario CA, Newark	508.9	127,916.8	12.8
3	Atlanta, Ontario CA, Newark	343.4	86,324.5	8.6
4	Atlanta, Chicago, Ontario CA, Newark	263.8	66,314.1	6.6
5	Atlanta, Chicago, Ontario CA, Denver, Newark	240.4	60,436.1	6.0

ggplot(pmedian, aes(n_dcs, wtd_miles)) +
  geom_line(color = "#1f77b4") + geom_point(size = 3, color = "#1f77b4") +
  labs(title = "Diminishing returns from adding DCs",
       x = "Number of DCs", y = "Demand-weighted miles to nearest DC") +
  theme_minimal(base_size = 13)

Discussion. Each DC has fixed costs (rent, labor, inventory). Where is the “knee” of the curve? How does the answer change if the emission factor or lead-time assumptions change?

8.3 Export data for Python, Tableau, and Sourcemap

Concept: CSV as a common language. Every tool in this lab can read a plain CSV file. Example: the same edges.csv becomes a network in R, a graph in Python, a flow map in Tableau, and a supply chain map in Sourcemap.

All three outside tools read the same CSV files.

dir.create("data", showWarnings = FALSE)

edges_geo <- edges %>%
  left_join(select(nodes, id, from_label = label, from_lon = lon, from_lat = lat), by = c("from" = "id")) %>%
  left_join(select(nodes, id, to_label = label,   to_lon = lon,   to_lat = lat),   by = c("to"   = "id"))

write_csv(nodes,                  "data/nodes.csv")
write_csv(edges_geo,              "data/edges.csv")
write_csv(national,               "data/national_weekly.csv")
write_csv(demand_weekly %>% select(id, label, week, demand), "data/customer_weekly_demand.csv")
write_csv(cust_avg,               "data/customer_avg_demand.csv")
write_csv(dc_reach,               "data/dc_reach.csv")
write_csv(stress,                 "data/stress_test.csv")

# A simple nodes/edges pair for Sourcemap (column names may need mapping on import)
write_csv(nodes %>% transmute(Name = paste(type, label, sep = " - "), Type = type,
                              Latitude = lat, Longitude = lon),
          "data/sourcemap_nodes.csv")
write_csv(edges_geo %>% transmute(From = paste(str_sub(from, 1, 1), from_label, sep = " - "),
                                  To = to, Tier = tier, Miles = miles),
          "data/sourcemap_edges.csv")

list.files("data")

## [1] "customer_avg_demand.csv"    "customer_weekly_demand.csv"
## [3] "dc_reach.csv"               "edges.csv"                 
## [5] "national_weekly.csv"        "nodes.csv"                 
## [7] "sourcemap_edges.csv"        "sourcemap_nodes.csv"       
## [9] "stress_test.csv"

9 Python: graph algorithms and a transportation model

Python is strong for optimization. Here we (a) recompute network metrics with networkx and (b) solve a transportation problem: ship each customer’s demand from the DCs at minimum total mile cost, with each DC limited to 28% of total demand.

Concept: shortest path. The route between two nodes with the smallest total distance. Example: if Houston to Dallas is 240 miles and Dallas to Atlanta is 780 miles, that route is 1,020 miles and wins if no other route is shorter.

Concept: transportation problem (linear programming). Decide how many units to ship from each source to each destination at minimum total cost, subject to capacity and demand constraints. Example: DCs A and B can each ship 100 units; customer X needs 80 and customer Y needs 90. A is cheaper for both, but its 100-unit limit forces the model to ship part of Y’s demand from B. That forced trade-off is exactly what “nearest DC” thinking misses.

First set the YAML parameter run_python: TRUE (top of this file) after installing the packages.

import pandas as pd, numpy as np, networkx as nx

nodes = pd.read_csv("data/nodes.csv")
edges = pd.read_csv("data/edges.csv")

G = nx.DiGraph()
for _, r in nodes.iterrows():
    G.add_node(r["id"], label=r["label"], type=r["type"])
for _, r in edges.iterrows():
    G.add_edge(r["from"], r["to"], miles=r["miles"])

# Betweenness centrality (weighted by miles)
bc = nx.betweenness_centrality(G, weight="miles")
top = sorted(bc.items(), key=lambda kv: kv[1], reverse=True)[:5]
print("Top 5 by betweenness:", [(G.nodes[n]["label"], round(v, 3)) for n, v in top])

# Shortest supplier-to-customer distance for every customer
sup  = nodes.loc[nodes.type == "Supplier", "id"]
cust = nodes.loc[nodes.type == "Customer", "id"]
rows = []
for c in cust:
    best = min((nx.shortest_path_length(G, s, c, weight="miles")
                for s in sup if nx.has_path(G, s, c)), default=np.nan)
    rows.append((G.nodes[c]["label"], round(best, 0)))
print(pd.DataFrame(rows, columns=["customer", "min_supplier_miles"]).head(10))

from math import radians, sin, cos, asin, sqrt
from scipy.optimize import linprog

def haversine(lon1, lat1, lon2, lat2):
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    a = sin((lat2 - lat1) / 2) ** 2 + cos(lat1) * cos(lat2) * sin((lon2 - lon1) / 2) ** 2
    return 3958.8 * 2 * asin(sqrt(a))          # Earth radius in miles

dcs = nodes[nodes.type == "DC"].reset_index(drop=True)
cu  = nodes[nodes.type == "Customer"].reset_index(drop=True)
dem = pd.read_csv("data/customer_avg_demand.csv").set_index("id").loc[cu.id, "avg_demand"].values

nd, nc = len(dcs), len(cu)
cost = np.array([[haversine(d.lon, d.lat, c.lon, c.lat) for _, c in cu.iterrows()]
                 for _, d in dcs.iterrows()])   # miles, shape (nd, nc)

cap  = np.full(nd, 0.28 * dem.sum())            # each DC can ship at most 28% of demand

# Decision variables x[i, j] flattened as i * nc + j
A_ub = np.kron(np.eye(nd), np.ones((1, nc)))    # DC capacity constraints
A_eq = np.kron(np.ones((1, nd)), np.eye(nc))    # every customer's demand is met

res = linprog(cost.flatten(), A_ub=A_ub, b_ub=cap, A_eq=A_eq, b_eq=dem,
              bounds=(0, None), method="highs")

flows = res.x.reshape(nd, nc)
print("Status:", res.message)
print("Optimal cost (unit-miles):", round(res.fun))
print("Nearest-DC cost (ignores capacity):", round((cost.min(axis=0) * dem).sum()))
print(pd.DataFrame({"DC": dcs.label, "units_shipped": flows.sum(axis=1).round(0),
                    "share": (flows.sum(axis=1) / dem.sum()).round(3)}))

Discussion. How does the optimal plan differ from “each customer buys from its nearest DC”? Which constraint is the reason?

10 Capstone assignment

Choose one product (e.g., coffee, solar panels, laptops, avocados) and:

Build a realistic 3-tier network (at least 3 suppliers, 2 plants, 3 DCs, 10 customer markets) with real coordinates.
Simulate or obtain a demand series and show one disruption (Demand and supply analysis).
Run a buffer analysis with at least two radii and identify underserved markets (Buffer and population-density analysis).
Compute centrality and run a stress test (Network analysis). Recommend one change that most improves resilience.
Build a Tableau dashboard and a Sourcemap map of the same network.
Publish the knitted HTML to RPubs and write a one-page executive summary: what is the biggest risk, and what would you do about it?

Grading component	Weight
Correct, reproducible R/Python analysis	35%
Tableau dashboard quality and interactivity	20%
Sourcemap map and tier-2 extension	15%
Insight and recommendations	20%
Clarity of writing and documentation	10%

11 References

Chopra, S., & Meindl, P. Supply Chain Management: Strategy, Planning, and Operation. Pearson.
Simchi-Levi, D., Kaminsky, P., & Simchi-Levi, E. Designing and Managing the Supply Chain. McGraw-Hill.
Hyndman, R. J., & Athanasopoulos, G. Forecasting: Principles and Practice. OTexts (free online).
Kolaczyk, E. D., & Csárdi, G. Statistical Analysis of Network Data with R. Springer.
Pebesma, E. (2018). Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal, 10(1), 439–446.
Csárdi, G., & Nepusz, T. (2006). The igraph software package for complex network research. InterJournal, Complex Systems, 1695.
Tableau documentation: spatial functions (MAKEPOINT, MAKELINE, BUFFER).
Sourcemap documentation (use the current help center for the import template).

Supply Chain Analytics Lab: Demand & Supply, Buffers, Networks, and Beyond

R + Python + Tableau + Sourcemap

Special Topics

November 10, 2026

1 Welcome and learning goals

1.1 Key vocabulary

1.2 Packages

2 Tableau: interactive maps and a decision dashboard

2.1 Build 1: Demand vs. supply (use `national_weekly.csv`)

2.2 Build 2: Network map (use `nodes.csv`)

2.3 Build 3: Flow map (use `edges.csv`)

2.4 Build 4: Interactive buffer map

2.5 Build 5: Dashboard

3 Sourcemap: multi-tier supply chain mapping

4 Our case: a simulated national supply chain

4.1 Spatial objects (`sf`)

5 Demand and supply analysis

5.1 Simulate three years of weekly demand

5.2 Simulate supply and a disruption

5.3 Inventory, shortfalls, and service level

5.4 Forecast the next quarter

6 Buffer and population-density analysis

6.1 Buffers around DCs

6.2 Reach of each DC

6.3 Map the buffers and find underserved customers

6.4 Optional: real population density with Census data

7 Network analysis

7.1 Build the network

7.2 Which nodes matter most?

7.3 Stress test: what if facilities fail?

8 Beyond the basics: inventory, network design, and emissions

8.1 Safety stock at each DC

8.2 Network design: how many DCs do we really need?

8.3 Export data for Python, Tableau, and Sourcemap

9 Python: graph algorithms and a transportation model

10 Capstone assignment

11 References

Supply Chain Analytics Lab: Demand & Supply, Buffers, Networks, and Beyond

R + Python + Tableau + Sourcemap

Special Topics

November 10, 2026

1 Welcome and learning goals

1.1 Key vocabulary

1.2 Packages

2 Tableau: interactive maps and a decision dashboard

2.1 Build 1: Demand vs. supply (use national_weekly.csv)

2.2 Build 2: Network map (use nodes.csv)

2.3 Build 3: Flow map (use edges.csv)

2.4 Build 4: Interactive buffer map

2.5 Build 5: Dashboard

3 Sourcemap: multi-tier supply chain mapping

4 Our case: a simulated national supply chain

4.1 Spatial objects (sf)

5 Demand and supply analysis

5.1 Simulate three years of weekly demand

5.2 Simulate supply and a disruption

5.3 Inventory, shortfalls, and service level

5.4 Forecast the next quarter

6 Buffer and population-density analysis

6.1 Buffers around DCs

6.2 Reach of each DC

6.3 Map the buffers and find underserved customers

6.4 Optional: real population density with Census data

7 Network analysis

7.1 Build the network

7.2 Which nodes matter most?

7.3 Stress test: what if facilities fail?

8 Beyond the basics: inventory, network design, and emissions

8.1 Safety stock at each DC

8.2 Network design: how many DCs do we really need?

8.3 Export data for Python, Tableau, and Sourcemap

9 Python: graph algorithms and a transportation model

10 Capstone assignment

11 References

2.1 Build 1: Demand vs. supply (use `national_weekly.csv`)

2.2 Build 2: Network map (use `nodes.csv`)

2.3 Build 3: Flow map (use `edges.csv`)

4.1 Spatial objects (`sf`)