Plant City Travel Demand Model

4-Step Model with Routed Skim & BPR Assignment

Author

Planning Team

Published

February 26, 2026

Scope — Internal Trips Only

This model currently includes internal-internal (I-I) trips only. External-internal (EI), internal-external (IE), and external-external (EE) trips are not yet incorporated. Assigned volumes therefore underrepresent traffic on corridors serving through and boundary movements.

# Set working directory to project root (QMD lives elsewhere)
knitr::opts_knit$set(root.dir = "D:/Renaissance Planning Group/Projects - TransportationNetworkStudy_25044/tech_group/data")

suppressPackageStartupMessages({
  library(tidyverse)
  library(sf)
  library(knitr)
  library(kableExtra)
  library(DT)
  library(plotly)
  library(tmap)
  library(leaflet)
  library(htmltools)
  library(cppRouting)
  library(igraph)
  library(flowmapblue)
  library(htmlwidgets)
  library(scales)
})

# Paths
zones_shp     <- "INTERIM/modeling/scenario_baseline_2020/zones.shp"
tripends_path <- "INTERIM/modeling/scenario_baseline_2020/trip_generation/trip_gen_internal_2020.csv"
tbrpm_path    <- "../code/shiny_apps/replica_flow_viewer/data/tbrpm_flows/flows_tbrpm.rds"
replica_path  <- "../code/shiny_apps/replica_flow_viewer/data/replica_flows/flows_replica.rds"
out_dir       <- "out_quick_gravity"

# Network prep outputs (from 00_prep_network.R — run that first!)
edges_rds      <- file.path(out_dir, "prepped_edges.rds")
connectors_rds <- file.path(out_dir, "prepped_connectors.rds")
nodes_rds      <- file.path(out_dir, "net_nodes.rds")
skim_csv       <- file.path(out_dir, "routed_skim.csv")
network_sf_rds <- file.path(out_dir, "network_sf.rds")

mapbox_token <- Sys.getenv("MAPBOX_API_TOKEN", unset = "")

purposes    <- c("HBW", "HBSH", "HBSR", "HBSC", "HBO", "NHBW", "NHBO")
hb_purposes <- c("HBW", "HBSH", "HBSR", "HBSC", "HBO")

# BPR parameters
bpr_alpha <- 0.15
bpr_beta  <- 4.0

Prerequisites

Run 00_prep_network.R before rendering this document. It builds the directed network, centroid connectors, and routed OD skim from the TBRPM network.

1 Introduction

This document implements a 4-step travel demand model for the Plant City study area (149 zones), following NCHRP 716/365 methodology and calibrated against TBRPM and Replica benchmark datasets.

Step	Method
Trip Generation	TBRPM-based rates, 7 purposes
Trip Distribution	NCHRP doubly-constrained gravity, gamma friction
Mode Split	Implicit (scale person-trips → TBRPM vehicle-trips)
Traffic Assignment	MSA with BPR volume-delay (\(\alpha=0.15\), \(\beta=4.0\))

The travel time skim is built by routing shortest paths on the actual network rather than using crow-fly distances. This replaced the original aggregated TBRPM skim which had an r=0.25 correlation with actual distances due to broken centroid connectors.

2 Step 0: Inputs

2.1 Zone System

Show code

zones_sf <- st_read(zones_shp, quiet = TRUE)
id_col   <- intersect(c("zone_id", "ZONE_ID", "TAZ", "taz", "ID", "id"),
                       names(zones_sf))[1]
zones_sf <- zones_sf %>% mutate(zone_id = as.integer(.data[[id_col]]))

zones_utm  <- zones_sf %>% st_transform(32617) %>%
  mutate(area_sqmi = as.numeric(st_area(.)) / 1609.34^2)
zones_4326 <- zones_sf %>% st_transform(4326)

Show code

pal_zone <- colorNumeric("Blues", zones_4326$zone_id)

leaflet(zones_4326) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor   = ~pal_zone(zone_id),
    fillOpacity = 0.25,
    color       = "#555555",
    weight      = 0.8,
    label       = ~paste0("Zone ", zone_id),
    highlightOptions = highlightOptions(weight = 2, color = "#333", bringToFront = TRUE)
  )

Figure 1: 149-zone system

Show code

zones_utm %>%
  st_drop_geometry() %>%
  summarise(`Zones` = n(),
            `Min Area (sq mi)` = round(min(area_sqmi), 2),
            `Median Area` = round(median(area_sqmi), 2),
            `Max Area` = round(max(area_sqmi), 2),
            `Total Area` = round(sum(area_sqmi), 1)) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 1: Zone system summary

Zones	Min Area (sq mi)	Median Area	Max Area	Total Area
149	0.03	0.93	35.27	297

2.2 Network

Show code

net_sf   <- readRDS(network_sf_rds)
edges    <- readRDS(edges_rds)
connectors <- readRDS(connectors_rds)
net_nodes  <- readRDS(nodes_rds)
edges_all  <- bind_rows(edges, connectors)

Show code

edges %>%
  filter(fac_class != "Snap Connector") %>%
  group_by(`Facility` = fac_class) %>%
  summarise(
    `Dir. Edges` = n(),
    `Total Miles` = round(sum(distance, na.rm = TRUE), 1),
    `Avg Speed (mph)` = round(mean(speed, na.rm = TRUE), 1),
    `Avg Hourly Cap` = format(round(mean(capacity, na.rm = TRUE)), big.mark = ","),
    .groups = "drop"
  ) %>%
  arrange(desc(`Dir. Edges`)) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 2: Network summary by facility class (excluding centroid connectors)

Facility	Dir. Edges	Total Miles	Avg Speed (mph)	Avg Hourly Cap
Collector	2144	1068.0	29.6	309
Divided Arterial	1038	433.8	38.4	457
Undivided Arterial	596	350.0	36.8	463
Toll Facility	142	134.7	44.8	1,633
Freeway	70	41.4	61.9	1,721
Ramp	58	22.0	29.3	1,350
One-Way Arterial	22	3.7	32.0	695

Show code

net_map <- net_sf %>%
  filter(fac_class != "Centroid Connector") %>%
  st_transform(4326)

fac_colors <- c(
  "Freeway"           = "#e41a1c",
  "Divided Arterial"  = "#ff7f00",
  "Undivided Arterial"= "#fdbf6f",
  "Collector"         = "#4daf4a",
  "Local"             = "#a6cee3",
  "Ramp"              = "#984ea3",
  "Snap Connector"    = "#999999"
)

fac_levels <- intersect(names(fac_colors), unique(net_map$fac_class))
pal_fac <- colorFactor(unname(fac_colors[fac_levels]), levels = fac_levels)

leaflet(net_map) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolylines(
    color   = ~pal_fac(fac_class),
    weight  = 1.5,
    opacity = 0.8,
    label   = ~paste0(fac_class, " | Speed: ", speed, " | Hourly Cap: ", capacity),
    highlightOptions = highlightOptions(weight = 3, opacity = 1)
  ) %>%
  addLegend("bottomright", pal = pal_fac, values = ~fac_class,
            title = "Facility Class", opacity = 0.9)

Figure 2: Road network colored by facility class

2.3 Routed Skim

Show code

skim <- read_csv(skim_csv, show_col_types = FALSE)

tripends <- read_csv(tripends_path, show_col_types = FALSE)
zones <- sort(intersect(unique(tripends$zone_id), unique(zones_sf$zone_id)))
zones <- zones[zones %in% unique(skim$origin)]
n_zones <- length(zones)

# Build travel time matrix
skim_filt <- skim %>% filter(origin %in% zones, destination %in% zones)
tt_mat <- matrix(NA_real_, n_zones, n_zones)
skim_idx <- skim_filt %>%
  mutate(i = match(origin, zones), j = match(destination, zones))
tt_mat[cbind(skim_idx$i, skim_idx$j)] <- skim_idx$travel_time
od_index <- expand_grid(origin = zones, destination = zones)

Show code

skim_filt %>%
  mutate(type = ifelse(origin == destination, "Intrazonal", "Interzonal")) %>%
  filter(is.finite(travel_time)) %>%
  group_by(Type = type) %>%
  summarise(Pairs = n(),
            `Min (min)` = round(min(travel_time), 1),
            `Median` = round(median(travel_time), 1),
            `Mean` = round(mean(travel_time), 1),
            `Max` = round(max(travel_time), 1),
            .groups = "drop") %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 3: Routed skim summary (free-flow network travel times)

Type	Pairs	Min (min)	Median	Mean	Max
Interzonal	22052	1.0	17.4	18	70.0
Intrazonal	149	0.4	2.5	3	15.4

Show code

inter_skim <- skim_filt %>%
  filter(origin != destination, is.finite(travel_time))

p <- ggplot(inter_skim, aes(x = travel_time)) +
  geom_histogram(binwidth = 2, fill = "#2c7fb8", color = "white", alpha = 0.8) +
  geom_vline(xintercept = mean(inter_skim$travel_time),
             color = "red", linetype = "dashed", linewidth = 0.8) +
  annotate("text", x = mean(inter_skim$travel_time) + 2,
           y = Inf, vjust = 2, hjust = 0, color = "red",
           label = sprintf("Mean: %.1f min", mean(inter_skim$travel_time))) +
  labs(x = "Travel Time (minutes)", y = "OD Pairs") +
  theme_minimal(base_size = 13)
ggplotly(p)

Figure 3: Distribution of interzonal routed travel times

3 Step 1: Trip Generation

Show code

tripends <- read_csv(tripends_path, show_col_types = FALSE)

Show code

tripgen_tbl <- map_dfr(purposes, function(p) {
  te <- tripends %>% filter(zone_id %in% zones)
  tibble(
    Purpose = p,
    Type = ifelse(p %in% hb_purposes, "HB", "NHB"),
    Productions = round(sum(te[[paste0("prod_", p)]], na.rm = TRUE)),
    Attractions = round(sum(te[[paste0("attr_", p)]], na.rm = TRUE)),
    `P/A Ratio` = round(Productions / Attractions, 3)
  )
})

tripgen_tbl %>%
  bind_rows(summarise(., Purpose = "**TOTAL**", Type = "—",
                      Productions = sum(Productions),
                      Attractions = sum(Attractions),
                      `P/A Ratio` = round(sum(Productions) / sum(Attractions), 3))) %>%
  kable(format = "html", format.args = list(big.mark = ","), escape = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(nrow(tripgen_tbl) + 1, bold = TRUE)

Table 4: Trip generation by purpose (daily person-trips)

Purpose	Type	Productions	Attractions	P/A Ratio
HBW	HB	90,725	76,274	1.189
HBSH	HB	37,447	36,345	1.030
HBSR	HB	104,833	116,663	0.899
HBSC	HB	58,665	49,845	1.177
HBO	HB	78,988	54,571	1.447
NHBW	NHB	18,698	18,973	0.986
NHBO	NHB	25,332	26,387	0.960
TOTAL	—	414,688	379,058	1.094

Show code

tripends %>%
  filter(zone_id %in% zones) %>%
  select(zone_id, starts_with("prod_"), starts_with("attr_")) %>%
  mutate(across(where(is.numeric), ~round(.))) %>%
  datatable(extensions = 'Buttons',
            options = list(pageLength = 15, scrollX = TRUE,
                           dom = 'Bfrtip', buttons = c('csv', 'excel')),
            rownames = FALSE)

Table 5: Trip ends by zone (download with buttons below)

Show code

prods <- tripends %>%
  filter(zone_id %in% zones) %>%
  mutate(total_prod = rowSums(across(starts_with("prod_")), na.rm = TRUE)) %>%
  select(zone_id, total_prod)

prods_sf <- zones_4326 %>%
  left_join(prods, by = "zone_id") %>%
  filter(!is.na(total_prod))

pal_prod <- colorQuantile("YlOrRd", prods_sf$total_prod, n = 6)

leaflet(prods_sf) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(
    fillColor   = ~pal_prod(total_prod),
    fillOpacity = 0.7,
    color       = "#444444",
    weight      = 0.5,
    label       = ~paste0("Zone ", zone_id, ": ", format(round(total_prod), big.mark = ","), " trips"),
    highlightOptions = highlightOptions(weight = 2, color = "#333", bringToFront = TRUE)
  ) %>%
  addLegend("bottomright", pal = pal_prod, values = ~total_prod,
            title = "Productions<br>(quantile)", opacity = 0.9)

Figure 4: Total trip productions by zone

4 Step 2: Trip Distribution

The gravity model uses a gamma friction function: \[F(t) = a \cdot t^b \cdot e^{ct}\] where \(c\) is calibrated per purpose to match TBRPM’s trip length distribution measured on the routed skim.

Show code

gravity_distribute <- function(P, A_target, F_mat, max_iter = 500, tol = 1e-4) {
  A <- A_target
  for (k in seq_len(max_iter)) {
    af    <- sweep(F_mat, 2, A, `*`)
    denom <- rowSums(af)
    denom[denom == 0] <- 1
    T_mat <- sweep(af, 1, P / denom, `*`)
    cs    <- colSums(T_mat)
    max_err <- max(abs(cs - A_target) / pmax(A_target, 1))
    if (max_err < tol) break
    A <- ifelse(cs > 0, A * A_target / cs, 0)
  }
  T_mat
}

pa_to_od <- function(PA_mat, is_hb) {
  if (is_hb) 0.5 * PA_mat + 0.5 * t(PA_mat) else PA_mat
}

build_friction_gamma <- function(tt_mat, a, b, cc) {
  tt <- pmax(tt_mat, 0.5)
  f_mat <- a * tt^b * exp(cc * tt)
  # Set Inf/NA cells to zero friction (unroutable pairs)
  f_mat[!is.finite(f_mat)] <- 0
  if (max(f_mat) > 0) f_mat <- f_mat / max(f_mat)
  f_mat
}

avg_trip_length <- function(OD_mat, tt_mat) {
  mask <- row(OD_mat) != col(OD_mat) & is.finite(tt_mat)
  sum(OD_mat[mask] * tt_mat[mask]) / sum(OD_mat[mask])
}

4.1 Calibration

Show code

# TBRPM benchmark
tbrpm <- readRDS(tbrpm_path) %>%
  filter(trip_type == "ALL") %>%
  transmute(origin, destination = dest, count)

tbrpm_total <- sum(tbrpm$count)

# Measure TBRPM avg trip length on routed skim
tbrpm_tt <- tbrpm %>%
  left_join(skim_filt, by = c("origin", "destination")) %>%
  filter(origin != destination, is.finite(travel_time))
tbrpm_avg_tt <- sum(tbrpm_tt$count * tbrpm_tt$travel_time) / sum(tbrpm_tt$count)

# Gamma base params (NCHRP 716 / TBRPM family)
gamma_a <- c(HBW = 28507, HBSH = 29013, HBSR = 63000,
             HBSC = 29013, HBO  = 29013, NHBW = 139173, NHBO = 219113)
gamma_b <- c(HBW = -0.020, HBSH = -0.400, HBSR = -0.700,
             HBSC = -0.400, HBO  = -0.350, NHBW = -0.800, NHBO = -0.900)

# Purpose targets — proportional to TBRPM avg, scaled by purpose
# These will be re-calibrated via bisection on c
purpose_target <- c(HBW = 11.0, HBSH = 7.5, HBSR = 8.5, HBSC = 6.0,
                    HBO = 8.0, NHBW = 8.5, NHBO = 7.0)

# Calibrate
max_iter <- 500; tol <- 1e-4

calibrated <- list()
cal_log <- list()

for (p in purposes) {
  pa <- tripends %>%
    filter(zone_id %in% zones) %>%
    arrange(match(zone_id, zones)) %>%
    transmute(P = replace_na(.data[[paste0("prod_", p)]], 0),
              A = replace_na(.data[[paste0("attr_", p)]], 0))

  P <- pa$P; A <- pa$A
  A_scaled <- A * sum(P) / sum(A)
  is_hb <- p %in% hb_purposes
  target <- purpose_target[p]

  c_lo <- -0.500; c_hi <- -0.005
  for (it in 1:30) {
    c_mid <- (c_lo + c_hi) / 2
    F_mat  <- build_friction_gamma(tt_mat, gamma_a[p], gamma_b[p], c_mid)
    PA_mat <- gravity_distribute(P, A_scaled, F_mat, max_iter, tol)
    OD_mat <- pa_to_od(PA_mat, is_hb)
    atl    <- avg_trip_length(OD_mat, tt_mat)
    if (atl > target) c_hi <- c_mid else c_lo <- c_mid
    if (abs(atl - target) < 0.1) break
  }

  calibrated[[p]] <- list(a = gamma_a[p], b = gamma_b[p], c = c_mid,
                          OD_mat = OD_mat, is_hb = is_hb)
  cal_log[[p]] <- tibble(
    Purpose = p, Type = ifelse(is_hb, "HB", "NHB"),
    `Target (min)` = target, `Achieved (min)` = round(atl, 1),
    a = gamma_a[p], b = gamma_b[p], c = round(c_mid, 4),
    Trips = round(sum(OD_mat)),
    `% Intrazonal` = round(sum(diag(OD_mat)) / sum(OD_mat) * 100, 1)
  )
}

Show code

bind_rows(cal_log) %>%
  kable(format = "html", format.args = list(big.mark = ",")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 6: Calibrated gamma friction parameters (routed skim)

Purpose	Type	Target (min)	Achieved (min)	a	b	c	Trips	% Intrazonal
HBW	HB	11.0	10.9	28,507	-0.02	-0.2216	90,725	9.5
HBSH	HB	7.5	8.7	29,013	-0.40	-0.5000	37,447	16.9
HBSR	HB	8.5	8.4	63,000	-0.70	-0.4381	104,833	20.7
HBSC	HB	6.0	9.1	29,013	-0.40	-0.5000	58,665	15.5
HBO	HB	8.0	8.1	29,013	-0.35	-0.4995	78,988	21.0
NHBW	NHB	8.5	8.5	139,173	-0.80	-0.2525	18,698	23.7
NHBO	NHB	7.0	7.1	219,113	-0.90	-0.5000	25,332	35.6

4.2 Scale to TBRPM & Build OD

Show code

raw_total <- sum(map_dbl(calibrated, ~sum(.x$OD_mat)))
scale_fac <- tbrpm_total / raw_total

od_all <- map_dfr(purposes, function(p) {
  OD <- calibrated[[p]]$OD_mat * scale_fac
  od_index %>% mutate(purpose = p,
                      travel_time = as.vector(t(tt_mat)),
                      trips = as.vector(t(OD)))
})

od_total <- od_all %>%
  summarise(trips = sum(trips), .by = c(origin, destination)) %>%
  left_join(od_index %>% mutate(travel_time = as.vector(t(tt_mat))),
            by = c("origin", "destination"))

write_csv(od_all,   file.path(out_dir, "od_by_purpose.csv"))
write_csv(od_total, file.path(out_dir, "od_all_purposes.csv"))

Show code

tibble(
  ` ` = c("Raw gravity (person-trips)", "TBRPM target (vehicle-trips)",
          "Scale factor", "TBRPM avg trip length on routed skim"),
  Value = c(format(round(raw_total), big.mark = ","),
            format(round(tbrpm_total), big.mark = ","),
            round(scale_fac, 4),
            paste0(round(tbrpm_avg_tt, 1), " min"))
) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)

Table 7: Trip scaling

	Value
Raw gravity (person-trips)	414,689
TBRPM target (vehicle-trips)	141,671
Scale factor	0.3416
TBRPM avg trip length on routed skim	9.7 min

4.3 Validation

Show code

replica <- readRDS(replica_path) %>%
  filter(trip_type == "Auto") %>%
  transmute(origin, destination = dest, replica = count)

compare <- od_total %>%
  transmute(origin, destination, gravity = trips) %>%
  full_join(tbrpm %>% rename(tbrpm = count), by = c("origin", "destination")) %>%
  full_join(replica, by = c("origin", "destination")) %>%
  replace_na(list(gravity = 0, tbrpm = 0, replica = 0)) %>%
  filter(origin != destination)

Show code

zone_cmp <- compare %>%
  group_by(origin) %>%
  summarise(across(c(gravity, tbrpm, replica), sum), .groups = "drop")

tibble(
  Comparison = c("Gravity vs TBRPM", "Gravity vs Replica", "TBRPM vs Replica"),
  `r (OD)` = round(c(cor(compare$gravity, compare$tbrpm),
                      cor(compare$gravity, compare$replica),
                      cor(compare$tbrpm, compare$replica)), 3),
  `r (Zone)` = round(c(cor(zone_cmp$gravity, zone_cmp$tbrpm),
                        cor(zone_cmp$gravity, zone_cmp$replica),
                        cor(zone_cmp$tbrpm, zone_cmp$replica)), 3)
) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 8: Correlation matrix (interzonal OD flows)

Comparison	r (OD)	r (Zone)
Gravity vs TBRPM	0.771	0.756
Gravity vs Replica	0.720	0.801
TBRPM vs Replica	0.779	0.858

Show code

compare_tt <- compare %>%
  left_join(skim_filt %>% filter(origin != destination),
            by = c("origin", "destination"))

tld <- bind_rows(
  compare_tt %>% transmute(Source = "Gravity", time = travel_time, trips = gravity),
  compare_tt %>% transmute(Source = "TBRPM", time = travel_time, trips = tbrpm),
  compare_tt %>% transmute(Source = "Replica Auto", time = travel_time, trips = replica)
) %>%
  filter(trips > 0, is.finite(time)) %>%
  mutate(bin = cut(time, breaks = c(0, 5, 10, 15, 20, 25, 30, 40, Inf),
                   include.lowest = TRUE)) %>%
  group_by(Source, bin) %>%
  summarise(trips = sum(trips), .groups = "drop") %>%
  group_by(Source) %>%
  mutate(share = trips / sum(trips) * 100)

p <- ggplot(tld, aes(x = bin, y = share, fill = Source)) +
  geom_col(position = position_dodge(0.8), alpha = 0.85) +
  scale_fill_manual(values = c("Gravity" = "#2c7fb8", "TBRPM" = "#d95f02",
                                "Replica Auto" = "#7570b3")) +
  labs(x = "Travel Time (min)", y = "Share (%)", fill = NULL) +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))
ggplotly(p)

Figure 5: Trip length distribution comparison (routed skim)

Show code

p <- ggplot(compare %>% filter(gravity > 0, tbrpm > 0),
            aes(x = tbrpm, y = gravity)) +
  geom_point(alpha = 0.2, size = 0.8, color = "#2c7fb8") +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  scale_x_log10(labels = comma) + scale_y_log10(labels = comma) +
  annotate("text", x = 1, y = max(compare$gravity) * 0.3, hjust = 0,
           label = sprintf("r = %.3f", cor(compare$gravity, compare$tbrpm)),
           color = "red", size = 5) +
  labs(x = "TBRPM (trips)", y = "Gravity (trips)") +
  theme_minimal(base_size = 13)
ggplotly(p)

Figure 6: Gravity vs TBRPM (interzonal, log scale)

4.4 Flowmap

Show code

locations <- zones_sf %>%
  filter(zone_id %in% zones) %>%
  st_transform(4326) %>% st_centroid() %>%
  mutate(lat = st_coordinates(.)[, 2], lon = st_coordinates(.)[, 1]) %>%
  st_drop_geometry() %>%
  transmute(id = as.character(zone_id), name = paste0("Zone ", zone_id), lat, lon)

flows <- od_total %>%
  filter(origin != destination, trips > 0) %>%
  transmute(origin = as.character(origin),
            dest = as.character(destination),
            count = round(trips))

flowmapblue(locations, flows,
            mapboxAccessToken = "pk.eyJ1IjoiamxlaG1hbjkzIiwiYSI6ImNrdmI2ZzQ2eWFrOXIycW1hZ3AwbGJnMWQifQ.HS4Uwi2cs5W5ZK0J7OsYkw",
            clustering = TRUE, darkMode = FALSE, animation = FALSE)

Figure 7: OD flow map (interactive, clustered)

5 Step 3: Mode Split

Implicit Mode Split

Scale factor of 0.3416 converts 414,689 person-trips to 141,671 vehicle-trips, capturing ~85% auto mode share and ~1.15 occupancy. A formal mode choice model is not included.

6 Step 4: Traffic Assignment

6.1 BPR Volume-Delay

\[t(V) = t_0 \left(1 + 0.15 \left(\frac{V}{C}\right)^4\right)\]

MSA iteration: \(V^{(n)} = V^{(n-1)} + \frac{1}{n}\left(V^{AON} - V^{(n-1)}\right)\)

Show code

bpr_time <- function(t0, vol, cap, a = 0.15, b = 4.0) {
  t0 * (1 + a * (pmax(vol, 0) / pmax(cap, 1))^b)
}

msa_iters <- 25
convergence_tol <- 0.001

# Prepare demand (scaled to TBRPM)
od_demand <- od_total %>%
  filter(origin != destination, trips > 0.5) %>%
  transmute(from = paste0("Z", origin),
            to   = paste0("Z", destination),
            demand = trips)

Show code

# CRITICAL: cppRouting deduplicates edges internally (keeps lowest cost per
# from→to pair). We must do the same so our edge table matches what cppRouting
# actually routes through. Otherwise path→volume lookup fails silently.

# Deduplicate: for each from→to pair, keep the edge with lowest free-flow time
edges_dedup <- edges_all %>%
  group_by(from, to) %>%
  slice_min(t0, n = 1, with_ties = FALSE) %>%
  ungroup() %>%
  mutate(eid = row_number())

# Now edge_key and the graph use the SAME edge set — no mismatches
edge_key <- edges_dedup %>%
  mutate(volume = 0, cost = t0)

# --- Daily capacity conversion ---
# Input capacities are hourly (HCM-based). Since the trip table is daily,
# convert to daily capacity using a peak-hour factor (K-factor).
# K ≈ 0.095 is typical for small urban areas (NCHRP 365 / FDOT).
K_factor <- 0.095
edge_key <- edge_key %>%
  mutate(hourly_cap = capacity,
         capacity   = capacity / K_factor)  # daily capacity ≈ hourly × 10.5

# Simple integer lookup: paste key → row index
edge_lu <- setNames(seq_len(nrow(edge_key)),
                     paste0(edge_key$from, "||", edge_key$to))

Show code

# Quick pre-assignment QC
n_real  <- sum(!edge_key$fac_class %in% c("Centroid Connector", "Snap Connector"))
n_cc    <- sum(edge_key$fac_class == "Centroid Connector")
n_total <- nrow(edge_key)
n_od    <- nrow(od_demand)
tot_dem <- sum(od_demand$demand)

# Test one path to verify loading works
test_graph <- makegraph(edge_key %>% transmute(from, to, dist = t0), directed = TRUE)
test_o <- od_demand$from[1]
test_d <- od_demand$to[1]
test_path <- get_path_pair(test_graph, from = test_o, to = test_d)[[1]]

# Verify every link in the test path is in our lookup
if (!is.null(test_path) && length(test_path) >= 2) {
  test_hits <- 0
  for (j in 1:(length(test_path) - 1)) {
    key <- paste0(test_path[j], "||", test_path[j + 1])
    if (!is.na(edge_lu[key])) test_hits <- test_hits + 1
  }
  path_match <- paste0(test_hits, "/", length(test_path) - 1, " edges matched")
} else {
  path_match <- "FAILED — no path found"
}

Pre-Assignment QC

Edges: 3,088 total (2,194 network + 894 connectors)
Demand: 10,844 OD pairs, 114,243 daily trips
Capacity: Hourly → daily (K = 0.095, factor ≈ 10.5×)
Path test (Z1 → Z2): 4/4 edges matched

Show code

origins <- unique(od_demand$from)
iter_log <- list()

for (iter in 1:msa_iters) {
  t0_iter <- Sys.time()

  # Build graph with current BPR costs (same edge set as edge_key)
  graph <- makegraph(
    edge_key %>% transmute(from, to, dist = cost),
    directed = TRUE
  )

  # All-or-nothing loading
  aon_vol <- numeric(nrow(edge_key))

  for (o in origins) {
    od_o <- od_demand %>% filter(from == o)
    if (nrow(od_o) == 0) next

    # get_path_pair needs equal-length from/to vectors
    paths <- tryCatch(
      get_path_pair(graph,
                    from = rep(o, nrow(od_o)),
                    to   = od_o$to),
      error = function(e) replicate(nrow(od_o), NULL, simplify = FALSE)
    )

    for (k in seq_along(paths)) {
      p <- paths[[k]]
      if (is.null(p) || length(p) < 2) next
      dem <- od_o$demand[k]
      for (j in 1:(length(p) - 1)) {
        idx <- edge_lu[paste0(p[j], "||", p[j + 1])]
        if (!is.na(idx)) aon_vol[idx] <- aon_vol[idx] + dem
      }
    }
  }

  # MSA blend
  lambda <- 1 / iter
  edge_key$volume <- edge_key$volume + lambda * (aon_vol - edge_key$volume)

  # BPR update
  edge_key$cost <- bpr_time(edge_key$t0, edge_key$volume, edge_key$capacity)

  # Relative gap
  aon_cost <- sum(aon_vol * edge_key$cost, na.rm = TRUE)
  cur_cost <- sum(edge_key$volume * edge_key$cost, na.rm = TRUE)
  rel_gap  <- ifelse(cur_cost > 0, abs(aon_cost - cur_cost) / cur_cost, 1)

  elapsed <- round(as.numeric(difftime(Sys.time(), t0_iter, units = "secs")), 1)
  loaded  <- sum(aon_vol > 0)   # links loaded in THIS iteration's AON
  max_vc  <- max(edge_key$volume / pmax(edge_key$capacity, 1))
  tot_aon <- round(sum(aon_vol))

  iter_log[[iter]] <- tibble(
    Iter = iter, `Rel Gap` = rel_gap,
    `AON Loaded Links` = loaded,
    `Cumul Loaded` = sum(edge_key$volume > 1),
    `Max V/C` = round(max_vc, 3),
    `AON Trips` = format(tot_aon, big.mark = ","),
    Seconds = elapsed
  )

  if (rel_gap < convergence_tol && iter > 3) break
}

iter_df <- bind_rows(iter_log)

Show code

# Post-assignment QC: how much demand actually got loaded?
total_aon_last <- sum(aon_vol)  # last iteration's AON
total_demand   <- sum(od_demand$demand)
pct_loaded     <- round(total_aon_last / total_demand * 100, 1)

# How many real (non-connector) edges carry volume?
real_with_vol <- edge_key %>%
  filter(!fac_class %in% c("Centroid Connector", "Snap Connector"),
         volume > 1) %>%
  nrow()

Assignment Loading Check

Demand: 114,243 trips
Last AON loaded: 1,188,084 trips (1040%)
Network links with volume > 1: 1,406
Iterations: 6

Show code

iter_df %>%
  mutate(`Rel Gap` = formatC(`Rel Gap`, format = "e", digits = 3)) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  scroll_box(height = "350px")

Table 9: MSA convergence log

Iter	Rel Gap	AON Loaded Links	Cumul Loaded	Max V/C	AON Trips	Seconds
1	0.000e+00	2121	2089	2.595	1,181,483	2.7
2	6.505e-03	2139	2115	1.298	1,226,831	2.4
3	2.749e-03	2121	2115	1.712	1,180,319	2.7
4	1.679e-03	2132	2115	1.845	1,183,323	2.6
5	1.153e-03	2133	2115	1.914	1,180,154	2.4
6	9.991e-04	2143	2115	1.871	1,188,084	2.9

Show code

p <- ggplot(iter_df, aes(x = Iter, y = `Rel Gap`)) +
  geom_line(color = "#d95f02", linewidth = 1) +
  geom_point(color = "#d95f02", size = 2) +
  geom_hline(yintercept = convergence_tol, linetype = "dashed", color = "gray50") +
  scale_y_log10() +
  labs(x = "Iteration", y = "Relative Gap (log)") +
  theme_minimal(base_size = 13)
ggplotly(p)

Figure 8: Convergence: relative gap by iteration

6.2 Results

Show code

# ---------- Identify real loaded edges (exclude connectors) ----------
real_loaded <- edge_key %>%
  filter(!fac_class %in% c("Centroid Connector", "Snap Connector"),
         volume > 0.5) %>%
  mutate(
    vc_ratio     = volume / pmax(capacity, 1),
    delay_factor = cost / pmax(t0, 0.01),
    speed_loaded = distance / pmax(cost / 60, 0.01)
  )

# ---------- Aggregate AB + BA to undirected links ----------
# KEY FIX: In 00_prep_network.R, transmute() evaluates sequentially, so
#   link_id = paste0(link_id, "_AB")  runs first,
#   orig_id = as.character(link_id)   picks up the MODIFIED value ("xxx_AB").
# net_sf$link_id has the ORIGINAL value ("xxx") with no suffix.
# Solution: strip the _AB / _BA suffix to recover the undirected link ID.

real_loaded <- real_loaded %>%
  mutate(undir_id = str_remove(orig_id, "_(AB|BA)$"))

link_volumes <- real_loaded %>%
  group_by(undir_id) %>%
  summarise(
    volume       = sum(volume),
    capacity     = sum(capacity),
    fac_class    = first(fac_class),
    distance     = first(distance),
    speed_ff     = first(speed),
    avg_cost     = weighted.mean(cost, w = pmax(volume, 0.001)),
    n_dirs       = n(),
    .groups      = "drop"
  ) %>%
  mutate(
    vc_ratio     = volume / pmax(capacity, 1),
    speed_loaded = distance / pmax(avg_cost / 60, 0.01),
    delay_factor = avg_cost / pmax(distance / pmax(speed_ff, 5) * 60, 0.01)
  )

# ---------- Geometry join ----------
link_geom <- net_sf %>%
  st_transform(4326) %>%
  mutate(link_id = as.character(link_id)) %>%
  select(link_id, geometry)

link_volumes <- link_volumes %>%
  mutate(undir_id = as.character(undir_id))

links_sf <- link_geom %>%
  inner_join(link_volumes, by = c("link_id" = "undir_id"))

# Diagnostic
if (nrow(links_sf) == 0) {
  warning("ZERO geometry matches — checking ID formats:")
  message("  net_sf link_id samples:  ", paste(head(link_geom$link_id, 5), collapse = ", "))
  message("  edge undir_id samples:   ", paste(head(link_volumes$undir_id, 5), collapse = ", "))
  message("  edge orig_id samples:    ", paste(head(real_loaded$orig_id, 5), collapse = ", "))
}

Show code

# How many links got geometry?
n_with_geom <- nrow(links_sf)
n_without   <- nrow(link_volumes) - n_with_geom

Geometry Join

1436 loaded links matched to network geometry (0 unmatched).

Show code

links_sf %>%
  st_drop_geometry() %>%
  group_by(Facility = fac_class) %>%
  summarise(Links = n(),
            `Total Vol` = format(round(sum(volume)), big.mark = ","),
            `Avg Vol` = round(mean(volume)),
            `Max Vol` = format(round(max(volume)), big.mark = ","),
            `Avg V/C` = round(mean(vc_ratio), 3),
            `Max V/C` = round(max(vc_ratio), 3),
            .groups = "drop") %>%
  arrange(desc(Links)) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 10: Volume summary by facility class

Facility	Links	Total Vol	Avg Vol	Max Vol	Avg V/C	Max V/C
Collector	755	370,874	491	4,175	0.160	1.871
Divided Arterial	288	311,877	1083	3,147	0.237	0.782
Undivided Arterial	215	106,465	495	1,903	0.125	0.482
Freeway	56	49,210	879	2,499	0.049	0.138
Toll Facility	53	65,587	1237	3,474	0.072	0.203
Ramp	47	21,712	462	1,531	0.034	0.153
One-Way Arterial	22	31,134	1415	2,169	0.195	0.311

Show code

links_sf %>%
  st_drop_geometry() %>%
  mutate(`V/C Range` = cut(vc_ratio,
    breaks = c(0, 0.5, 0.75, 0.9, 1.0, 1.25, Inf),
    labels = c("<0.50", "0.50–0.75", "0.75–0.90",
               "0.90–1.00", "1.00–1.25", ">1.25"),
    include.lowest = TRUE)) %>%
  count(`V/C Range`, .drop = FALSE) %>%
  mutate(`%` = round(n / sum(n) * 100, 1)) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)

Table 11: V/C distribution

V/C Range	n	%
<0.50	1359	94.6
0.50–0.75	63	4.4
0.75–0.90	8	0.6
0.90–1.00	0	0.0
1.00–1.25	2	0.1
>1.25	4	0.3

Show code

vmt <- sum(links_sf$volume * links_sf$distance, na.rm = TRUE)
vht <- sum(links_sf$volume * links_sf$avg_cost / 60, na.rm = TRUE)

tibble(
  Measure = c("Daily VMT", "Daily VHT", "Network Avg Speed", "Loaded Links",
              "Total Daily Trips"),
  Value = c(paste(format(round(vmt), big.mark = ","), "mi"),
            paste(format(round(vht), big.mark = ","), "hrs"),
            paste(round(vmt / pmax(vht, 0.01), 1), "mph"),
            format(nrow(links_sf), big.mark = ","),
            format(round(tbrpm_total), big.mark = ","))
) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped"), full_width = FALSE)

Table 12: System performance

Measure	Value
Daily VMT	334,730 mi
Daily VHT	9,687 hrs
Network Avg Speed	34.6 mph
Loaded Links	1,436
Total Daily Trips	141,671

6.3 Maps

Show code

# Background: all real network links (no connectors) for context
net_bg <- net_sf %>%
  filter(!fac_class %in% c("Centroid Connector", "Snap Connector")) %>%
  st_transform(4326)

# Round values for clean popups / labels
links_sf <- links_sf %>%
  mutate(
    vol_round   = round(volume),
    cap_round   = round(capacity),
    vc_round    = round(vc_ratio, 2),
    speed_ld_rd = round(speed_loaded, 1),
    speed_ff_rd = round(speed_ff, 1),
    # Popup HTML for all maps
    popup_html = paste0(
      "<strong>", fac_class, "</strong><br>",
      "Daily Volume: ", format(vol_round, big.mark = ","), "<br>",
      "Daily Capacity: ", format(cap_round, big.mark = ","), "<br>",
      "V/C: ", vc_round, "<br>",
      "FF Speed: ", speed_ff_rd, " mph<br>",
      "Loaded Speed: ", speed_ld_rd, " mph"
    )
  )

# Line weight scaling helper: map volume to 1–8 px range
vol_range <- range(links_sf$vol_round, na.rm = TRUE)
links_sf$lwd_vol <- scales::rescale(links_sf$vol_round, to = c(1.5, 8), from = vol_range)

Show code

pal_vol <- colorNumeric("YlOrRd", links_sf$vol_round)

leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolylines(
    data    = net_bg,
    color   = "#d0d0d0",
    weight  = 0.6,
    opacity = 0.5,
    group   = "Network"
  ) %>%
  addPolylines(
    data    = links_sf,
    color   = ~pal_vol(vol_round),
    weight  = ~lwd_vol,
    opacity = 0.85,
    popup   = ~popup_html,
    label   = ~paste0(fac_class, ": ", format(vol_round, big.mark = ",")),
    highlightOptions = highlightOptions(weight = 6, opacity = 1, bringToFront = TRUE),
    group   = "Volumes"
  ) %>%
  addLegend(
    "bottomright", pal = pal_vol, values = links_sf$vol_round,
    title = "Daily Volume", opacity = 0.9
  ) %>%
  addLayersControl(
    overlayGroups = c("Network", "Volumes"),
    options = layersControlOptions(collapsed = FALSE)
  )

Figure 9: Assigned daily volumes

Show code

# LOS-style categorical V/C
links_sf <- links_sf %>%
  mutate(vc_cat = cut(vc_ratio,
    breaks = c(0, 0.5, 0.75, 0.9, 1.0, 1.25, Inf),
    labels = c("<0.50", "0.50-0.75", "0.75-0.90",
               "0.90-1.00", "1.00-1.25", ">1.25"),
    include.lowest = TRUE))

vc_colors <- c(
  "<0.50"     = "#2166ac",
  "0.50-0.75" = "#67a9cf",
  "0.75-0.90" = "#fdd49e",
  "0.90-1.00" = "#ef8a62",
  "1.00-1.25" = "#b2182b",
  ">1.25"     = "#67001f"
)

vc_levels <- levels(links_sf$vc_cat)
pal_vc <- colorFactor(unname(vc_colors[vc_levels]), levels = vc_levels, ordered = TRUE)

leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolylines(
    data    = net_bg,
    color   = "#d0d0d0",
    weight  = 0.6,
    opacity = 0.5,
    group   = "Network"
  ) %>%
  addPolylines(
    data    = links_sf,
    color   = ~pal_vc(vc_cat),
    weight  = 2.5,
    opacity = 0.85,
    popup   = ~popup_html,
    label   = ~paste0(fac_class, " V/C: ", vc_round),
    highlightOptions = highlightOptions(weight = 5, opacity = 1, bringToFront = TRUE),
    group   = "V/C"
  ) %>%
  addLegend(
    "bottomright", pal = pal_vc, values = links_sf$vc_cat,
    title = "V/C Ratio", opacity = 0.9
  ) %>%
  addLayersControl(
    overlayGroups = c("Network", "V/C"),
    options = layersControlOptions(collapsed = FALSE)
  )

Figure 10: Volume-to-capacity ratio

Show code

# Fixed speed breaks (mph)
links_sf <- links_sf %>%
  mutate(spd_cat = cut(speed_ld_rd,
    breaks = c(0, 10, 20, 30, 40, 50, Inf),
    labels = c("<10", "10-20", "20-30", "30-40", "40-50", ">50"),
    include.lowest = TRUE))

spd_colors <- c(
  "<10"   = "#d73027",
  "10-20" = "#fc8d59",
  "20-30" = "#fee08b",
  "30-40" = "#d9ef8b",
  "40-50" = "#91cf60",
  ">50"   = "#1a9850"
)

spd_levels <- levels(links_sf$spd_cat)
pal_spd <- colorFactor(unname(spd_colors[spd_levels]), levels = spd_levels, ordered = TRUE)

leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolylines(
    data    = net_bg,
    color   = "#d0d0d0",
    weight  = 0.6,
    opacity = 0.5,
    group   = "Network"
  ) %>%
  addPolylines(
    data    = links_sf,
    color   = ~pal_spd(spd_cat),
    weight  = 2.5,
    opacity = 0.85,
    popup   = ~popup_html,
    label   = ~paste0(fac_class, ": ", speed_ld_rd, " mph"),
    highlightOptions = highlightOptions(weight = 5, opacity = 1, bringToFront = TRUE),
    group   = "Speed"
  ) %>%
  addLegend(
    "bottomright", pal = pal_spd, values = links_sf$spd_cat,
    title = "Loaded Speed (mph)", opacity = 0.9
  ) %>%
  addLayersControl(
    overlayGroups = c("Network", "Speed"),
    options = layersControlOptions(collapsed = FALSE)
  )

Figure 11: Congested speed

6.4 Top Loaded Corridors

Show code

links_sf %>%
  st_drop_geometry() %>%
  arrange(desc(volume)) %>%
  head(25) %>%
  transmute(
    `Link ID`      = link_id,
    Facility       = fac_class,
    Miles          = round(distance, 2),
    `FF Speed`     = speed_ff_rd,
    Volume         = vol_round,
    `Daily Cap`    = cap_round,
    `V/C`          = vc_round,
    `Loaded Speed` = speed_ld_rd
  ) %>%
  kable(format = "html", format.args = list(big.mark = ",")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 13: Top 25 loaded links

Link ID	Facility	Miles	FF Speed	Volume	Daily Cap	V/C	Loaded Speed
20616_20582	Collector	0.14	25	4,175	2,232	1.87	8.8
20582_20616	Collector	0.14	25	4,094	2,232	1.83	9.3
20589_20616	Collector	0.14	25	3,835	2,232	1.72	10.8
20616_20589	Collector	0.14	25	3,818	2,232	1.71	10.9
20713_20730	Toll Facility	0.14	44	3,474	17,137	0.20	13.7
20729_20713	Toll Facility	0.14	44	3,474	17,137	0.20	13.7
20578_20577	Collector	0.15	30	3,439	2,968	1.16	14.9
8774_7224	Collector	0.21	32	3,302	3,895	0.85	20.7
7224_8774	Collector	0.21	32	3,219	3,895	0.83	20.7
7273_8163	Divided Arterial	0.17	38	3,147	4,026	0.78	17.0
8163_7273	Divided Arterial	0.17	38	3,141	4,026	0.78	17.0
20581_20578	Divided Arterial	0.14	38	3,131	4,384	0.71	14.3
20582_20581	Divided Arterial	0.14	38	3,131	4,384	0.71	14.3
20577_20578	Collector	0.15	30	3,018	2,968	1.02	14.9
7315_8742	Divided Arterial	0.31	38	2,906	4,384	0.66	30.8
8742_7273	Divided Arterial	0.19	38	2,906	4,384	0.66	19.4
7273_8742	Divided Arterial	0.19	38	2,887	4,384	0.66	19.4
8742_7315	Divided Arterial	0.31	38	2,887	4,384	0.66	30.8
20578_20581	Divided Arterial	0.14	38	2,687	4,384	0.61	14.3
20581_20582	Divided Arterial	0.14	38	2,687	4,384	0.61	14.3
20746_20747	Collector	0.18	32	2,668	3,989	0.67	17.6
8140_8193	Divided Arterial	0.50	38	2,570	4,384	0.59	37.3
8193_8225	Divided Arterial	0.25	38	2,570	4,384	0.59	25.4
8225_8249	Divided Arterial	0.25	38	2,570	4,384	0.59	24.8
20573_20746	Collector	0.18	32	2,568	3,989	0.64	17.6

6.5 Link-Level Data

Show code

links_sf %>%
  st_drop_geometry() %>%
  transmute(link_id, fac_class,
            miles    = round(distance, 2),
            ff_speed = speed_ff_rd,
            volume   = vol_round,
            capacity = cap_round,
            vc       = vc_round,
            loaded_speed = speed_ld_rd,
            delay    = round(delay_factor, 3)) %>%
  arrange(desc(volume)) %>%
  datatable(extensions = 'Buttons',
            options = list(pageLength = 20, scrollX = TRUE,
                           dom = 'Bfrtip', buttons = c('csv', 'excel')),
            rownames = FALSE) %>%
  formatRound(columns = c("volume", "capacity"), digits = 0) %>%
  formatRound(columns = c("vc", "delay"), digits = 3)

Table 14: All loaded links (sortable, downloadable)

Show code

# Save assigned links CSV
write_csv(
  links_sf %>% st_drop_geometry(),
  file.path(out_dir, "assigned_links.csv")
)

7 Summary

Show code

tibble(
  Step = c("Trip Generation", "Trip Distribution",
           "Mode Split", "Traffic Assignment"),
  Method = c("TBRPM rates, 7 purposes",
             "Doubly-constrained gravity, gamma friction, routed skim",
             paste0("Implicit (factor = ", round(scale_fac, 4), ")"),
             paste0("MSA × ", nrow(iter_df), " iters, BPR (α=0.15, β=4)")),
  Result = c(
    paste0(format(round(raw_total), big.mark = ","), " person-trips"),
    paste0("r = ", round(cor(compare$gravity, compare$tbrpm), 3), " vs TBRPM"),
    paste0(format(round(tbrpm_total), big.mark = ","), " vehicle-trips"),
    paste0("VMT = ", format(round(vmt), big.mark = ","),
           " | Speed = ", round(vmt / pmax(vht, 0.01), 1), " mph")
  )
) %>%
  kable(format = "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Table 15: 4-Step model results

Step	Method	Result
Trip Generation	TBRPM rates, 7 purposes	414,689 person-trips
Trip Distribution	Doubly-constrained gravity, gamma friction, routed skim	r = 0.771 vs TBRPM
Mode Split	Implicit (factor = 0.3416)	141,671 vehicle-trips
Traffic Assignment	MSA × 6 iters, BPR (α=0.15, β=4)	VMT = 334,730 \| Speed = 34.6 mph

Output files in out_quick_gravity/:

File	Description
`od_by_purpose.csv`	OD matrix by purpose (7 purposes, scaled)
`od_all_purposes.csv`	Total OD matrix
`assigned_links.csv`	Link volumes, V/C, loaded speeds
`routed_skim.csv`	Network-routed free-flow travel times
`prepped_edges.rds`	Directed edge table