Blood Transshipment modelling
Aiden Scholey
2026-09-17
Intro and loading data
Blood is redistributed and shipped across Canada. Some cities operate in a surplus of certain blood types and others are constantly in deficits. This sparked the thought of creating an LP model to attempt to identify the best possible outcome for blood shipping in Canada. The issue was that I couldn’t find any Canadian data online- so I decided to go with data from Australia - and examine them instead.
The dataset I initially found on Kaggle (https://www.kaggle.com/datasets/ahmadmehmod/blood-donor-dataset) contained synthetic blood donation records for 10,000 blood donation events across Australia. The wrinkle was that they were divided evenly across 8 major cities in Australia, and also evenly split across the 8 blood groups (12.5% for each city and blood group). This is not representative of the general population in Australia. So the first task at hand was re balancing the dataset so that it included more O and A blood donors (who are the most common blood types based on genetic occurrence).
This rebalanced dataset was then loaded into R.
setwd("C:/Users/ascholey/Desktop/R Files/Blood Donor")
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'tibble' was built under R version 4.4.3## Warning: package 'tidyr' was built under R version 4.4.3## Warning: package 'readr' was built under R version 4.4.3## Warning: package 'purrr' was built under R version 4.4.3## Warning: package 'dplyr' was built under R version 4.4.3## Warning: package 'stringr' was built under R version 4.4.3## Warning: package 'forcats' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(lpSolve)## Warning: package 'lpSolve' was built under R version 4.4.2donors2 = read.csv("australian_blood_donor_dataset_rebalanced.csv")Continuing the rebalancing
Then we continue rebalancing the data by reshuffling the donor records to be better representative of the population distribution in Australia.
In 2a we introduce the population estimates for Australia’s eight capital cities as reference points.
In 2b we introduce the population_share table which represents the estimated prevalence of each ABO/RhD blood group in the Australian population.
In 2b.ii we make an important distinction by introducing donor pool share. I didnt want to assume that the donor population in Australia is the same as the general population distribution of blood types. So I found data from a publication (Vox Sanguinis 2024) of blood donor prevalence in Australia. This is criticial, because donor recruiting patterns are not necessarily representative of population genetics.
In 2c we alter the donor records in the dataset (donors2) and they decide where they live/donate. This was done by calculating city population/total AU capital city population and then multiplied that by the number of synthetic donors (and then dealt with rounding). Lastly we shuffle the donor-city assignments before being attached to individual donor records.
In 2d we add a unique element to this analysis/dataset, replacing the population-share-equals-demand-share assumption. It factors in an “emergency O- usage” parameter algebraically so that modeled O- demand lands on the real observed value (~16% of hospital orders). This treats O- as disproportionately used in emergency/unknown type situations not just routine group specific matching.
This was calculated as: (1-E)(pop’n_share) + E(emergency_oneg_share) = observed_demand, and solved for E.
# 2a. Real capital city populations, 30 June 2024 (ABS)
city_population <- tribble(
~city, ~population,
"Sydney", 5557200,
"Melbourne", 5350700,
"Brisbane", 2780100,
"Adelaide", 1469200,
"Perth", 2384400,
"Hobart", 254900,
"Darwin", 152500,
"Canberra", 473900
)
# 2b. National ABO/RhD population prevalence (Vox Sanguinis 2024 National Prevelance)
population_share <- tribble(
~blood_group, ~pop_share,
"O+", 0.380,
"O-", 0.070,
"A+", 0.320,
"A-", 0.060,
"B+", 0.120,
"B-", 0.020,
"AB+", 0.040,
"AB-", 0.010
)
#2b.ii AU Donor pool composition (Vox Sanguinis 2024 n = 566k)
# reflects actual donor recruitment patterns rather than pop'n genetics
donor_pool_share <- tribble(
~blood_group, ~donor_share,
"O+", 0.387,
"O-", 0.107,
"A+", 0.297,
"A-", 0.060,
"B+", 0.098,
"B-", 0.017,
"AB+", 0.029,
"AB-", 0.005
)
# 2c. City allocation
n_donors <- nrow(donors2)
# Calculate target number of donors per city
city_allocation <- city_population %>%
mutate(
population_share = population / sum(population),
target_donors_raw = population_share * n_donors,
target_donors = floor(target_donors_raw)
)
# Distribute remaining donors according to largest fractional remainder
remaining_donors <- n_donors - sum(city_allocation$target_donors)
city_allocation <- city_allocation %>%
arrange(desc(target_donors_raw - target_donors)) %>%
mutate(
target_donors = target_donors +
if_else(row_number() <= remaining_donors, 1L, 0L)
)
# Create city assignment
set.seed(123)
city_assignment <- rep(
city_allocation$city,
city_allocation$target_donors
)
# Shuffle city assignments
city_assignment <- sample(city_assignment)
# Shuffle donor records independently
donors2 <- donors2 %>%
slice_sample(prop = 1)
# Assign cities
donors2$city <- city_assignment
# CHECK
donors2 %>%
count(city) %>%
left_join(
city_allocation %>%
select(city, population, population_share, target_donors),
by = "city"
)# 2d. Modelled RBC demand share, incorporating emergency group-O usage
recipient_share <- population_share %>%
rename(recipient_share = pop_share)
# Emergency-release fractions based on 2024 South Australia MHP O usage (n=700)
emerg_oneg_share <- 0.45
emerg_opos_share <- 0.55
# Calibrate the emergency group-O fraction so overall O- demand matches
# the observed Australian hospital-order share (~16%)
o_neg_population_share <- population_share %>%
filter(blood_group == "O-") %>%
pull(pop_share)
observed_o_neg_demand <- 0.16
emergency_group_o_share <-
(observed_o_neg_demand - o_neg_population_share) /
(emerg_oneg_share - o_neg_population_share)
demand_share <- recipient_share %>%
mutate(
routine_demand = recipient_share * (1 - emergency_group_o_share),
emergency_demand = case_when(
blood_group == "O-" ~ emergency_group_o_share * emerg_oneg_share,
blood_group == "O+" ~ emergency_group_o_share * emerg_opos_share,
TRUE ~ 0
),
demand_share = routine_demand + emergency_demand
) %>%
select(blood_group, demand_share)Resampling blood groups and donor pool population
In 3a we resample the blood_group, by overwriting the datasets blood-type mix by drawing from a new variable donor_pool_share instead of the population_share. This is where the donor-recruitiment skew data gets built into the simulated data.
In 3b a chi-square test is run to confirm the re sampling worked as planned, observed counts in donors2 stastiscaly match the target donor_pool_share distribution p = ~0.5 meaning no meaningful deviation.
Lastly, we close this section by printing a plot of the donor pool vs population prevalence comparison. This visualizes the evidence of real world recruitment skew (O- donor pool ~1.5x its population share).
# Is the synthetic donor pool's blood-type mix realistic?
# 3a. Resample blood_group to match real donor pop'n
set.seed(123)
n_donors = nrow(donors2)
donors2 = donors2 %>%
mutate(
blood_group = sample(
donor_pool_share$blood_group,
size = n_donors,
replace = TRUE,
prob = donor_pool_share$donor_share
)
)
# 3b. Is the synthetic donor pool's blood-type mix realistic?
observed_counts <- donors2 %>% count(blood_group, name = "n") %>% arrange(blood_group)
expected_props <- donor_pool_share %>%
arrange(blood_group) %>%
pull(donor_share)
expected_counts <- expected_props / sum(expected_props) * sum(observed_counts$n)
chisq_result <- chisq.test(x = observed_counts$n, p = expected_props / sum(expected_props))
print(chisq_result)##
## Chi-squared test for given probabilities
##
## data: observed_counts$n
## X-squared = 6.2994, df = 7, p-value = 0.5053comparison_df <- donor_pool_share %>%
left_join(population_share, by = "blood_group") %>%
mutate(
skew_ratio = donor_share / pop_share # >1 = over-recruited relative to population
)
print(comparison_df)## # A tibble: 8 × 4
## blood_group donor_share pop_share skew_ratio
## <chr> <dbl> <dbl> <dbl>
## 1 O+ 0.387 0.38 1.02
## 2 O- 0.107 0.07 1.53
## 3 A+ 0.297 0.32 0.928
## 4 A- 0.06 0.06 1
## 5 B+ 0.098 0.12 0.817
## 6 B- 0.017 0.02 0.85
## 7 AB+ 0.029 0.04 0.725
## 8 AB- 0.005 0.01 0.5ggplot(
comparison_df %>%
pivot_longer(c(donor_share, pop_share), names_to = "series", values_to = "share"),
aes(x = blood_group, y = share, fill = series)
) +
geom_col(position = "dodge") +
scale_y_continuous(labels = scales::percent) +
scale_fill_manual(
values = c(donor_share = "#c0392b", pop_share = "#2c3e50"),
labels = c(donor_share = "Donor pool (Vox Sanguinis 2024)", pop_share = "Community prevalence (Vox Sanguinis 2024)")
) +
labs(
title = "Australian donor pool vs. community prevalence by blood group",
subtitle = "Donor pool skews toward O- and O+ relative to population genetics",
x = "Blood group", y = "Share", fill = NULL
) +
theme_minimal()
City x Blood-Type Gaps
In section 4 the concept of fair share is defined. This represents what each city-bloodtype cell should hold if national supply were distributed purely by population share x demand share. Comparing that to the actual supply illustrates geographic mismatches. The heatmap here is visual evidence displaying that supply and demand dont line up by city.
# 4. City x blood-type supply, and each city's "fair share" of demand
supply <- donors2 %>%
group_by(city, blood_group) %>%
summarise(pints_supplied = sum(pints_donated), .groups = "drop")
total_supply_national <- sum(supply$pints_supplied)
total_population <- sum(city_population$population)
fair_share <- city_population %>%
crossing(demand_share) %>%
mutate(
city_pop_share = population / total_population,
fair_share_pints = city_pop_share * demand_share * total_supply_national
) %>%
select(city, blood_group, fair_share_pints)
gap <- supply %>%
full_join(fair_share, by = c("city", "blood_group")) %>%
mutate(
pints_supplied = replace_na(pints_supplied, 0),
gap = pints_supplied - fair_share_pints # positive = surplus, negative = deficit
)
print(gap %>% arrange(gap) %>% head(10)) # biggest deficits## # A tibble: 10 × 5
## city blood_group pints_supplied fair_share_pints gap
## <chr> <chr> <int> <dbl> <dbl>
## 1 Melbourne O- 10577 17479. -6902.
## 2 Sydney O- 12190 18154. -5964.
## 3 Sydney O+ 43436 47684. -4248.
## 4 Melbourne O+ 42680 45912. -3232.
## 5 Perth O- 4851 7789. -2938.
## 6 Brisbane O- 6446 9082. -2636.
## 7 Brisbane O+ 21846 23855. -2009.
## 8 Perth O+ 18864 20459. -1595.
## 9 Adelaide O+ 11483 12606. -1123.
## 10 Adelaide O- 3872 4799. -927.print(gap %>% arrange(desc(gap)) %>% head(10)) # biggest surpluses## # A tibble: 10 × 5
## city blood_group pints_supplied fair_share_pints gap
## <chr> <chr> <int> <dbl> <dbl>
## 1 Melbourne A+ 34775 26679. 8096.
## 2 Sydney A+ 33970 27708. 6262.
## 3 Brisbane A+ 17557 13862. 3695.
## 4 Perth A+ 14726 11889. 2837.
## 5 Sydney A- 7732 5195. 2537.
## 6 Adelaide B+ 3625 2747. 878.
## 7 Adelaide A+ 8137 7325. 812.
## 8 Melbourne A- 5751 5002. 749.
## 9 Brisbane B+ 5944 5198. 746.
## 10 Melbourne B+ 10708 10005. 703.ggplot(gap, aes(x = blood_group, y = city, fill = gap)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "#c0392b", mid = "white", high = "#27ae60", midpoint = 0) +
labs(
title = "Supply gap vs. demand-weighted fair share, by city x blood group",
subtitle = "Red = deficit relative to fair share | Green = surplus",
x = "Blood group", y = NULL, fill = "Pints\n(gap)"
) +
theme_minimal()
Reallocation LP (by blood type)
Section 5 is an initial attempt at a model. We have 8 independent transportation problems (1 x each ABORh), moving surplus-city pints to deficit-city pints, minimizing distance. This LP does not allow substitution between compatible blood types. This limitation is addressed in the later full LP.
# 5. Reallocation optimization: move surplus pints to cover deficits
cities <- c("Sydney", "Melbourne", "Brisbane", "Adelaide", "Perth", "Hobart", "Darwin", "Canberra")
dist_km <- matrix(
c(
0, 880, 920, 1400, 3900, 1050, 3150, 290, # Sydney
880, 0, 1670, 730, 3400, 600, 3140, 470, # Melbourne
920, 1670, 0, 1600, 3600, 1800, 2850, 970, # Brisbane
1400, 730, 1600, 0, 2130, 1160, 2620, 1050, # Adelaide
3900, 3400, 3600, 2130, 0, 3600, 2650, 3000, # Perth
1050, 600, 1800, 1160, 3600, 0, 3800, 900, # Hobart
3150, 3140, 2850, 2620, 2650, 3800, 0, 3140, # Darwin
290, 470, 970, 1050, 3000, 900, 3140, 0 # Canberra
),
nrow = 8, byrow = TRUE, dimnames = list(cities, cities)
)
allocate_one_blood_group <- function(bg, gap_df, dist_mat) {
cell <- gap_df %>% filter(blood_group == bg)
surplus <- cell %>% filter(gap > 0) %>% select(city, amount = gap)
deficit <- cell %>% filter(gap < 0) %>% mutate(amount = -gap) %>% select(city, amount)
if (nrow(surplus) == 0 || nrow(deficit) == 0) return(NULL)
cost_mat <- dist_mat[surplus$city, deficit$city, drop = FALSE]
# Balance the transportation problem (lpSolve's transport solver needs
# equal total supply and demand) by capping at the smaller total.
total_move <- min(sum(surplus$amount), sum(deficit$amount))
supply_vec <- surplus$amount * (total_move / sum(surplus$amount))
demand_vec <- deficit$amount * (total_move / sum(deficit$amount))
solution <- lp.transport(
cost.mat = cost_mat,
row.signs = rep("<=", nrow(cost_mat)),
row.rhs = supply_vec,
col.signs = rep("<=", ncol(cost_mat)),
col.rhs = demand_vec
)
moves <- as.data.frame(solution$solution)
colnames(moves) <- deficit$city
moves$from_city <- surplus$city
moves <- moves %>%
pivot_longer(-from_city, names_to = "to_city", values_to = "pints_moved") %>%
filter(pints_moved > 0.5) %>%
mutate(blood_group = bg, distance_km = map2_dbl(from_city, to_city, ~dist_mat[.x, .y]))
moves
}
reallocation_plan <- map_dfr(unique(gap$blood_group), allocate_one_blood_group, gap_df = gap, dist_mat = dist_km)
print(reallocation_plan %>% arrange(desc(pints_moved)))## # A tibble: 0 × 5
## # ℹ 5 variables: from_city <chr>, to_city <chr>, pints_moved <dbl>,
## # blood_group <chr>, distance_km <dbl>reallocation_plan %>%
summarise(
total_pints_moved = sum(pints_moved),
total_pint_km = sum(pints_moved * distance_km)
) %>%
print()## # A tibble: 1 × 2
## total_pints_moved total_pint_km
## <dbl> <dbl>
## 1 0 0Demand vs. Population comparison
In Section 6, demand_share vs raw population share was revisited with a repeat of the city-gap analysis using a finalized demand_share. I also troubleshoot for N/A cells (cities with small populations and zero donors of rare types like AB-)
#6. Population vs Modelled Demand
demand_comparison <- recipient_share %>%
left_join(demand_share, by = "blood_group") %>%
mutate(
demand_vs_population = demand_share / recipient_share
)
print(demand_comparison)## # A tibble: 8 × 4
## blood_group recipient_share demand_share demand_vs_population
## <chr> <dbl> <dbl> <dbl>
## 1 O+ 0.38 0.420 1.11
## 2 O- 0.07 0.16 2.29
## 3 A+ 0.32 0.244 0.763
## 4 A- 0.06 0.0458 0.763
## 5 B+ 0.12 0.0916 0.763
## 6 B- 0.02 0.0153 0.763
## 7 AB+ 0.04 0.0305 0.763
## 8 AB- 0.01 0.00763 0.763# and plot the results
ggplot(
demand_comparison %>%
select(blood_group, recipient_share, demand_share) %>%
pivot_longer(
c(recipient_share, demand_share),
names_to = "series",
values_to = "share"
),
aes(x = blood_group, y = share, fill = series)
) +
geom_col(position = "dodge") +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Australian population vs. modelled RBC demand",
subtitle = "Model incorporates emergency group-O utilization",
x = "Blood group",
y = "Share",
fill = NULL
) +
theme_minimal()
## Calculate supply by group and city
supply = donors2 %>%
group_by(city, blood_group) %>%
summarise(
pints_supplied = sum(pints_donated),
.groups = "drop"
)
total_supply_national = sum(supply$pints_supplied)
total_population = sum(city_population$population)
fair_share = city_population %>%
crossing(demand_share) %>%
mutate(
city_popn_share = population/total_population,
fair_share_pints =
city_popn_share *
demand_share *
total_supply_national
) %>%
select(city, blood_group, fair_share_pints)
## Then calculate the ratio between supply and demand
supply_gap <- supply %>%
full_join(
fair_share,
by = c("city", "blood_group")
) %>%
mutate(
pints_supplied = replace_na(pints_supplied, 0),
gap = pints_supplied - fair_share_pints,
supply_ratio = pints_supplied / fair_share_pints
)
supply_gap <- supply_gap %>%
mutate(
log_supply_ratio = log2(supply_ratio)
)
## Then plot
ggplot(
supply_gap,
aes(
x = blood_group,
y = city,
fill = log_supply_ratio
)
) +
geom_tile(color = "white") +
scale_fill_gradient2(
low = "#b2182b",
mid = "white",
high = "#1a9850",
midpoint = 0
) +
labs(
title = "Synthetic blood supply relative to fair share",
subtitle = "White = fair share | Red = deficit | Green = surplus",
x = "Blood group",
y = NULL,
fill = "log₂(Supply /\nFair share)"
) +
theme_minimal()
##### trouble shooting the code for the plot
city_check <- donors2 %>%
count(city, name = "donors") %>%
left_join(
city_population %>%
mutate(
population_share = population / sum(population)
),
by = "city"
) %>%
mutate(
donor_share = donors / sum(donors),
donor_vs_population = donor_share / population_share
)
### Checking the plot for Na's because Hobart and Darwin have N/A for AB- in the plot
supply_gap %>%
filter(blood_group == "AB-", city %in% c("Hobart", "Darwin"))## Now run a national blood group summary as discussed
national_blood_balance <- supply_gap %>%
group_by(blood_group) %>%
summarise(
supply = sum(pints_supplied),
fair_share = sum(fair_share_pints),
supply_ratio = supply / fair_share,
gap = supply - fair_share,
.groups = "drop"
) %>%
arrange(supply_ratio)
national_blood_balanceCompatibility Matrix
In section 7 we build a compatibility matrix based on standard RBC donor-recipient compatibility rules. This will later become the backbone for substitution logic in the final LP.
#7. Building a compatibility logic matrix for transfusion
rbc_compatibility <- tribble(
~donor, ~recipient, ~compatible,
# O-
"O-", "O-", TRUE,
"O-", "O+", TRUE,
"O-", "A-", TRUE,
"O-", "A+", TRUE,
"O-", "B-", TRUE,
"O-", "B+", TRUE,
"O-", "AB-", TRUE,
"O-", "AB+", TRUE,
# O+
"O+", "O-", FALSE,
"O+", "O+", TRUE,
"O+", "A-", FALSE,
"O+", "A+", TRUE,
"O+", "B-", FALSE,
"O+", "B+", TRUE,
"O+", "AB-", FALSE,
"O+", "AB+", TRUE,
# A-
"A-", "O-", FALSE,
"A-", "O+", FALSE,
"A-", "A-", TRUE,
"A-", "A+", TRUE,
"A-", "B-", FALSE,
"A-", "B+", FALSE,
"A-", "AB-", TRUE,
"A-", "AB+", TRUE,
# A+
"A+", "O-", FALSE,
"A+", "O+", FALSE,
"A+", "A-", FALSE,
"A+", "A+", TRUE,
"A+", "B-", FALSE,
"A+", "B+", FALSE,
"A+", "AB-", FALSE,
"A+", "AB+", TRUE,
# B-
"B-", "O-", FALSE,
"B-", "O+", FALSE,
"B-", "A-", FALSE,
"B-", "A+", FALSE,
"B-", "B-", TRUE,
"B-", "B+", TRUE,
"B-", "AB-", TRUE,
"B-", "AB+", TRUE,
# B+
"B+", "O-", FALSE,
"B+", "O+", FALSE,
"B+", "A-", FALSE,
"B+", "A+", FALSE,
"B+", "B-", FALSE,
"B+", "B+", TRUE,
"B+", "AB-", FALSE,
"B+", "AB+", TRUE,
# AB-
"AB-", "O-", FALSE,
"AB-", "O+", FALSE,
"AB-", "A-", FALSE,
"AB-", "A+", FALSE,
"AB-", "B-", FALSE,
"AB-", "B+", FALSE,
"AB-", "AB-", TRUE,
"AB-", "AB+", TRUE,
# AB+
"AB+", "O-", FALSE,
"AB+", "O+", FALSE,
"AB+", "A-", FALSE,
"AB+", "A+", FALSE,
"AB+", "B-", FALSE,
"AB+", "B+", FALSE,
"AB+", "AB-", FALSE,
"AB+", "AB+", TRUE
)
compatibility_matrix = rbc_compatibility %>%
select(donor, recipient, compatible) %>%
pivot_wider(
names_from = recipient,
values_from = compatible
)
compatibility_matrixrbc_compatibility %>%
filter(compatible) %>%
group_by(recipient) %>%
summarise(
compatible_donors = paste(donor, collapse = ", "),
.groups = "drop"
)National Compatibility Analysis
In section 8 I attempted to ask “how much compatible blood exists for each recipient type nationally”. This analysis was flawed by the fact that I double counted O- and O+ supply across every recipient type that compatible with them, since a single unit cant cover multiple recipients. Additionally, O- reservation logic was defined here, carving out enough O- supply for O- recipients first, so other types only draw from true leftover surplus.
#8. National Compatibility Analysis
# Here we calculate total # of pints avail for each donor group nationally
national_supply = donors2 %>%
group_by(blood_group) %>%
summarise(
supply = sum(pints_donated),
.groups = "drop"
)
# Recall that we have national demand as demand_share from section 2
# In the model as of now total demand is equal to available supply
total_demand = sum(national_supply$supply)
# Convert demand % from section 2 to actual pint quantity
national_demand = demand_share %>%
mutate(
demand_pints = demand_share * total_demand
)
# Then create every possible donor recipient combo
compatible_pairs = rbc_compatibility %>%
filter(compatible == TRUE)
# Now we need to determine the total compatible supply available to each recipient
# consists of identifying all compatible blood groups and adding national inventory together
compatible_supply = compatible_pairs %>%
left_join(
national_supply,
by = c("donor" = "blood_group")
) %>%
group_by(recipient) %>%
summarise(
compatible_supply = sum(supply),
.groups = "drop"
)
# Add the modelled demand for each recipient group and calculate the
# amount of compatible inventory that could theoretically cover it.
compatibility_summary <- national_demand %>%
select(blood_group, demand_pints) %>%
left_join(
compatible_supply,
by = c("blood_group" = "recipient")
) %>%
mutate(
potential_coverage = compatible_supply / demand_pints
)
compatibility_summary#attempting to protect o- recipients first here
#should pull down the inflated covereage ratios of A+, B+, AB+ since lost access to o-
o_neg_supply_total <- national_supply %>% filter(blood_group == "O-") %>% pull(supply)
o_neg_demand_pints <- national_demand %>% filter(blood_group == "O-") %>% pull(demand_pints)
o_neg_surplus_for_sharing <- max(o_neg_supply_total - o_neg_demand_pints, 0)
# Everyone else only sees leftover O- surplus, not the whole O- pool
national_supply_reserved <- national_supply %>%
mutate(supply = if_else(blood_group == "O-", o_neg_surplus_for_sharing, supply))
compatible_supply_reserved <- compatible_pairs %>%
left_join(national_supply_reserved, by = c("donor" = "blood_group")) %>%
group_by(recipient) %>%
summarise(compatible_supply = sum(supply), .groups = "drop")
compatibility_summary_reserved <- national_demand %>%
select(blood_group, demand_pints) %>%
left_join(compatible_supply_reserved, by = c("blood_group" = "recipient")) %>%
mutate(
# O- recipients get their full reserved pool back, not the shared leftover
compatible_supply = if_else(blood_group == "O-", o_neg_supply_total, compatible_supply),
potential_coverage = compatible_supply / demand_pints
)
compatibility_summary_reservedcomparison <- compatibility_summary %>%
select(blood_group, unreserved_coverage = potential_coverage) %>%
left_join(
compatibility_summary_reserved %>%
select(blood_group, reserved_coverage = potential_coverage),
by = "blood_group"
)
comparisonThe Final LP - Sections A-B
Now we begin to create the LP. In the first 2 sections (A,B) we define supply/demand tables, and enumarate variables
# -----------------------------------------------------------------------
# A. Supply and demand tables (unchanged from v1)
# -----------------------------------------------------------------------
lp_supply <- donors2 %>%
group_by(city, blood_group) %>%
summarise(supply = sum(pints_donated), .groups = "drop") %>%
rename(donor_type = blood_group)
lp_demand <- fair_share %>%
rename(recipient_type = blood_group, demand = fair_share_pints) %>%
mutate(demand_row_id = row_number())
# -----------------------------------------------------------------------
# B. Shipment variables (unchanged from v1)
# -----------------------------------------------------------------------
compatible_pairs <- rbc_compatibility %>% filter(compatible)
variable_index <- compatible_pairs %>%
rename(donor_type = donor, recipient_type = recipient) %>%
select(donor_type, recipient_type) %>%
crossing(origin = cities, dest = cities) %>%
mutate(
distance_km = map2_dbl(origin, dest, ~ dist_km[.x, .y]),
var_id = row_number()
)
n_ship_vars <- nrow(variable_index)The Final LP - Sections C-D
This section documents the infeasibility: hard demand constraints made the LP unsolvable initially because some blood types are structurally short nationally (B-, A-). This is fixed here by creating one shortfall variable for each demand cell. This variable is essentially saying “How many units of this particular city x recipient group demand cannot be fulfilled?”. In section D we define the objective, which is to minimize transportation cost, plus minimize shortfall penalty. O- was shortfall penalized 5x more than other types, an attempt to prioritize the idea of protect O- first as a cost weight rather than a rigid rule.
# C. NEW: one shortfall (slack) variable per demand cell, shortfall[i] = how much demand cell i falls short of fair_share
n_slack_vars <- nrow(lp_demand)
n_vars_total <- n_ship_vars + n_slack_vars
SHORTFALL_PENALTY <- 50000
o_neg_penalty_multiplier <- 5
# D. Objective: distance cost for shipments + penalty cost for shortfalls
slack_penalties <- if_else(
lp_demand$recipient_type == "O-",
SHORTFALL_PENALTY * o_neg_penalty_multiplier,
SHORTFALL_PENALTY
)
objective <- c(variable_index$distance_km, slack_penalties)The Final LP - Sections E-G
In this section of the LP, constraints are defined such as supply constraints capping what each city can ship out, demand constraints with slack to target what each city needs, the O- reservation constraint caps how much O- can go to non-O- recipients nationally, based on true surplus after O- demand is met.
# E. Supply constraints
supply_keys <- lp_supply %>% mutate(row_id = row_number())
supply_constraints <- matrix(0, nrow = nrow(supply_keys), ncol = n_vars_total)
for (i in seq_len(nrow(supply_keys))) {
matches <- variable_index$donor_type == supply_keys$donor_type[i] &
variable_index$origin == supply_keys$city[i]
supply_constraints[i, seq_len(n_ship_vars)][matches] <- 1
}
supply_rhs <- supply_keys$supply
supply_signs <- rep("<=", nrow(supply_keys))
# F. Demand constraints - NOW AN EQUALITY WITH SLACK:
# shipped_in + shortfall == demand
# (equivalently: shipped_in >= demand - shortfall, with shortfall >= 0
# handled automatically since lp() defaults to non-negative variables)-
demand_constraints <- matrix(0, nrow = nrow(lp_demand), ncol = n_vars_total)
for (i in seq_len(nrow(lp_demand))) {
ship_matches <- variable_index$recipient_type == lp_demand$recipient_type[i] &
variable_index$dest == lp_demand$city[i]
demand_constraints[i, seq_len(n_ship_vars)][ship_matches] <- 1
# This row's own slack variable gets coefficient 1
demand_constraints[i, n_ship_vars + lp_demand$demand_row_id[i]] <- 1
}
demand_rhs <- lp_demand$demand
demand_signs <- rep("=", nrow(lp_demand))
# G. O- hard reservation
o_neg_supply_total <- lp_supply %>% filter(donor_type == "O-") %>% pull(supply) %>% sum()
o_neg_demand_total <- lp_demand %>% filter(recipient_type == "O-") %>% pull(demand) %>% sum()
o_neg_surplus_cap <- max(o_neg_supply_total - o_neg_demand_total, 0)
o_neg_to_others <- variable_index$donor_type == "O-" & variable_index$recipient_type != "O-"
o_neg_constraint_row <- c(as.numeric(o_neg_to_others), rep(0, n_slack_vars))The Final LP - Section H-I
In the last section we solve the LP and unpack the results. lp_result$status == 0 confirming the feasibility. shipment_plan and shortfall_report are 2 key outputs. One showing the optimal routes and one showing where demand could not be fully met.
# H. Assemble and solve
const_mat <- rbind(supply_constraints, demand_constraints, o_neg_constraint_row)
const_dir <- c(supply_signs, demand_signs, "<=")
const_rhs <- c(supply_rhs, demand_rhs, o_neg_surplus_cap)
lp_result <- lp(
direction = "min",
objective.in = objective,
const.mat = const_mat,
const.dir = const_dir,
const.rhs = const_rhs
)
lp_result$status # should now be 0## [1] 0lp_result$objval## [1] 5672087342# I. Unpack results: shipment plan AND unmet-demand report
all_solution <- lp_result$solution
shipment_plan <- variable_index %>%
mutate(pints_shipped = all_solution[seq_len(n_ship_vars)]) %>%
filter(pints_shipped > 0.5) %>%
arrange(desc(pints_shipped))
shortfall_report <- lp_demand %>%
mutate(unmet_pints = all_solution[n_ship_vars + demand_row_id]) %>%
filter(unmet_pints > 0.5) %>%
arrange(desc(unmet_pints))
print(shipment_plan)## # A tibble: 84 × 7
## donor_type recipient_type origin dest distance_km var_id pints_shipped
## <chr> <chr> <chr> <chr> <dbl> <int> <dbl>
## 1 O+ O+ Sydney Sydney 0 1728 43436
## 2 O+ O+ Melbourne Melbour… 0 1710 42680
## 3 A+ A+ Sydney Sydney 0 320 27708.
## 4 A+ A+ Melbourne Melbour… 0 302 26679.
## 5 O+ O+ Brisbane Brisbane 0 1674 21846
## 6 O+ O+ Perth Perth 0 1719 18864
## 7 A+ A+ Brisbane Brisbane 0 266 13862.
## 8 O- O- Sydney Sydney 0 1408 12190
## 9 A+ A+ Perth Perth 0 311 11889.
## 10 O+ O+ Adelaide Adelaide 0 1665 11483
## # ℹ 74 more rowsprint(shortfall_report)## # A tibble: 17 × 5
## city recipient_type demand demand_row_id unmet_pints
## <chr> <chr> <dbl> <int> <dbl>
## 1 Melbourne O- 17479. 47 6902.
## 2 Sydney O- 18154. 63 5964.
## 3 Sydney O+ 47684. 64 4176.
## 4 Melbourne O+ 45912. 48 3232.
## 5 Perth O- 7789. 55 2938.
## 6 Brisbane O- 9082. 15 2636.
## 7 Brisbane O+ 23855. 16 2009.
## 8 Perth O+ 20459. 56 1595.
## 9 Adelaide O- 4799. 7 927.
## 10 Adelaide O+ 12606. 8 790.
## 11 Canberra O- 1548. 23 605.
## 12 Melbourne B- 1667. 45 289.
## 13 Hobart O- 833. 39 148.
## 14 Hobart O+ 2187. 40 120.
## 15 Darwin O- 498. 31 113.
## 16 Perth B- 743. 53 18.0
## 17 Darwin B- 47.5 29 13.5shipment_plan %>%
summarise(total_pints_shipped = sum(pints_shipped), total_pint_km = sum(pints_shipped * distance_km))shortfall_report %>%
group_by(recipient_type) %>%
summarise(total_unmet = sum(unmet_pints), .groups = "drop") %>%
arrange(desc(total_unmet))shortfall_report = shortfall_report %>%
mutate(pct_unmet = unmet_pints/demand)
shortfall_reportshortfall_vs_population <- shortfall_report %>%
left_join(city_population, by = "city")
ggplot(shortfall_vs_population, aes(x = population, y = pct_unmet, color = recipient_type)) +
geom_point(size = 3) +
geom_text(aes(label = city), vjust = -0.8, size = 3, show.legend = FALSE) +
scale_x_continuous(labels = scales::comma) +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Shortfall severity vs. city population",
subtitle = "Testing whether smaller cities face proportionally worse shortages",
x = "City population",
y = "% of demand unmet",
color = "Recipient type"
) +
theme_minimal()
Results & Discussion
The scatter plot above tests a natural hypothesis: that smaller cities, with less local donor infrastructure, would face proportionally worse blood shortages after optimal reallocation. The data doesn’t support this, Canberra (population ~474k) and Melbourne (population ~5.35M) show nearly identical O- shortfall rates (~39%), while Hobart and Darwin, the two smallest cities, actually perform comparatively well.
Blood type, not city size, is the dominant driver of shortfall severity. O- consistently shows 17-40% unmet demand across every city regardless of population, while O+ consistently sits at 5-9% in the same cities. Tracing this further back to each city’s local O- supply-to-demand ratio (calculated independently, before any LP reallocation) explains the pattern much more directly than population size does.
Limitations
This data is purely synthetic and was mapped to be similar to Australia’s national population. Real Life Blood data would make this more realistic and representative of the constraints on the Australian blood system. I also didnt incorporate any measure of time/blood shelf-life into this analysis. Future LP’s could attempt to build in the limited shelf-life of RBC’s which further constrains the blood system. Lastly, I made arbitrary shortfall penalty weights. One could argue that these could be changed and they arent based on actual values.
Conclusions
This model is a first step towards creating something that represents my home country (Canada) and our own blood transport system constraints. I was only able to find Australian blood donor data initially (at least in a large scale data file), so finding Canadian data (or making it) and then tuning this model to fit our domestic stats would be a logical next step. Thanks for reading!