54  Inventory Analytics and Optimisation

Author

Bongo Adi

Note📋 Learning Objectives
  • Understand inventory costs and the inventory-service level tradeoff
  • Apply ABC-XYZ segmentation to prioritise inventory management
  • Calculate economic order quantity (EOQ) and optimal reorder points
  • Incorporate demand and lead time uncertainty into safety stock models
  • Formulate and solve linear programming for supply chain networks
  • Simulate stochastic inventory policies and compare performance
  • Implement inventory optimisation with Nigerian pharmaceutical case study

54.1 The Invisible Drain on Business Profits

Somewhere in a warehouse in Aba, right now, there are 40,000 boxes of a nutritional supplement that nobody wants. The supplier delivered them six months ago because the purchasing manager — worried about running out — ordered three months’ worth of buffer stock. But demand was lower than expected, a competitor launched a similar product, and now the boxes sit and wait. Each day that passes, those boxes consume warehouse space (which could hold faster-moving products), tie up working capital (money borrowed at 20% annual interest that is not earning a return), and inch closer to their expiry date. In three months, 15,000 of those boxes will be beyond their sell-by date and will be destroyed. The total cost: ₦12 million written off, plus months of carrying costs.

Now imagine the opposite scenario in the same city. A pharmacy discovers that a critical antibiotic — the only one effective against a strain of bacteria circulating in the region — is out of stock. Patients who need it today cannot get it. Some will go to a competitor and never return. Others will suffer harm from delayed treatment. The pharmacy loses the sale, loses the customer, and loses its reputation for reliability. The cost of that stockout: incalculable in human terms, and hundreds of thousands of naira in lost revenue and customer lifetime value.

Between these two extremes — too much and too little — lies the fundamental challenge of inventory management. Every business that holds physical goods must solve this puzzle: how much to stock, when to reorder, how much safety buffer to keep against uncertainty. Get it right and the business runs efficiently with happy customers and healthy cash flow. Get it wrong in either direction and the cost is real: either cash locked in dusty shelves or customers lost to stockouts.

Analytics turns this from guesswork into a science. In this chapter, we build the quantitative tools that let a business find the optimal inventory level: models that balance the cost of holding stock against the cost of running out, that account for uncertainty in demand and lead time, and that prioritise scarce management attention across hundreds or thousands of products.

54.1.1 The Cost of Inventory: Why It Is Never Free

It is tempting to think that holding a large amount of stock is “safe” — if you have plenty, you can never run out. But inventory is never free. Every unit of unsold stock costs money, and those costs add up to a significant fraction of a company’s total expenses.

Holding (carrying) costs include:

  • Capital cost: The money tied up in inventory is not available for other uses. If the company’s cost of capital is 20% per year (typical in Nigeria given high lending rates), and inventory is worth ₦100 million, the opportunity cost is ₦20 million per year — the return the business would have earned by investing that money elsewhere.
  • Storage cost: Warehouses charge rent, consume electricity, and require staff. A 5,000 sq ft warehouse in Lagos might cost ₦500,000 per month. If 60% of that space is inventory, the monthly storage cost is ₦300,000.
  • Obsolescence and expiry: Medicines, food, and cosmetics expire. Electronics become outdated. Fashion goes out of style. For FMCG (Fast-Moving Consumer Goods) companies, expiration loss runs at 2–5% of inventory value annually — a constant drain.
  • Insurance and shrinkage: Goods in storage need to be insured against fire and theft. And despite best efforts, some inventory disappears through theft, damage, or spoilage — adding another 2–3% annually.

Total holding cost as a percentage of inventory value typically runs 20–35% per year in Nigerian conditions — higher than in many markets because of elevated borrowing costs, less reliable storage infrastructure, and greater shrinkage rates.

Ordering costs (the cost of placing and receiving a new order):

  • Fixed cost per order (paperwork, supplier negotiation, receiving inspection): ₦50,000–₦200,000 per order
  • Variable cost per unit: minimal, but sometimes exists for bulk handling or freight

Stockout costs (the cost of running out):

  • Lost sales: the margin on units you could have sold but did not. For a ₦5,000 product with 30% margin, each lost sale costs ₦1,500.
  • Expedited shipping: if you can emergency-source the product, air freight vs. sea freight might cost ₦5,000–₦50,000 extra per unit.
  • Lost customers: a customer who finds an empty shelf may never return.
  • Reputational damage: impossible to fully quantify but very real — especially for pharmacies, hospitals, and businesses where reliability is a core promise.

In Nigeria specifically, carrying costs are elevated for several reasons beyond high interest rates: security costs for guarding warehouses, foreign exchange risk (imported raw materials locked in inventory lose value if the naira weakens), and port delays that mean goods arriving from overseas are in transit — and therefore unavailable — for unpredictably long periods. A Nigerian pharmaceutical distributor carrying ₦500 million in inventory faces ₦100–150 million in annual carrying costs, creating a powerful incentive to reduce stock while maintaining service levels.

Ordering costs (for procurement/replenishment): - Fixed cost per order (administrative, shipping): ₦50k–₦200k per order - Variable cost (per unit): ₦0–₦10 depending on product

Stockout costs (when inventory runs out): - Lost sales (if customer buys elsewhere): margin on lost unit × lost quantity - Expedited shipping (air freight vs sea): ₦5,000–₦50,000 premium per unit - Goodwill loss & customer churn: harder to quantify but significant

In Nigerian context, carrying costs are elevated: security costs (guarding warehouses), foreign exchange risk (imported raw materials lose value if naira weakens), port delays (products tied up longer in customs), and theft (shrinkage in northern Nigeria estimated at 5–10% annually). A Nigerian pharmaceutical distributor carrying ₦500m inventory faces ₦100–150m annual carrying costs, justifying aggressive inventory reduction.

Note📘 Theory: Inventory Cost Tradeoffs

Holding cost increases with inventory level; ordering cost decreases (fewer orders). Total cost = holding + ordering cost. The Economic Order Quantity (EOQ) minimises total cost by balancing these tradeoffs. At EOQ, marginal holding cost = marginal ordering cost.

Tip🔑 Key Formula: Total Inventory Cost

\[\text{Total Cost} = \text{Holding Cost} + \text{Ordering Cost}\] \[= \frac{Q}{2} \times h + \frac{D}{Q} \times s\]

where \(Q\) = order quantity, \(h\) = holding cost per unit-year, \(D\) = annual demand, \(s\) = fixed cost per order.

Minimising: \(\text{EOQ} = \sqrt{\frac{2Ds}{h}}\)

Caution📝 Section 49.1 Review Questions

  1. A Nigerian pharmaceutical distributor carries ₦500m inventory. Holding costs are 25% annually, ordering cost per order is ₦100k. What are the three main types of holding costs? How does each scale with inventory level?

  2. Explain why holding costs are higher in Nigeria than in developed markets. Give three specific reasons related to supply chain environment.

  3. Total inventory cost = ₦(Q/2)h + ₦(D/Q)s. If holding cost doubles from ₦50 to ₦100 per unit-year, does total cost increase or decrease at the current order quantity? Why?

  4. A stockout of a high-margin product costs ₦500 in lost profit; a stockout of a low-margin product costs ₦50. How should the company’s target service levels differ for these two products?

  5. Define the three inventory costs (holding, ordering, stockout) and rank them by magnitude for: (a) a fast-moving consumer good (high volume, low margin), (b) a specialty medicine (low volume, high margin).

54.2 ABC-XYZ Analysis: Prioritise Your Inventory

Not all SKUs (stock keeping units) deserve equal attention. A pharmaceutical distributor carries 500 SKUs but might find that 80 SKUs represent 80% of inventory value. Managing expensive, slow-moving items differently from cheap, fast-moving items optimises working capital.

ABC Segmentation (by value): - A items (~10–20% of SKUs, 70–80% of value): High-value, often slow-moving. Expensive drugs, specialty medicines. - B items (~20–30% of SKUs, 15% of value): Mid-value. - C items (~50–70% of SKUs, 5–10% of value): Low-value, high-volume. Paracetamol tablets, bandages.

XYZ Segmentation (by demand variability): - X items: Stable demand (coefficient of variation CV < 0.5). Predictable, easy to forecast. - Y items: Variable demand (CV 0.5–1.0). Seasonal or promotional swings. - Z items: Erratic demand (CV > 1.0). Sporadic purchases, hard to forecast.

Combined ABC-XYZ Matrix (3×3):

X (Stable) Y (Variable) Z (Erratic)
A (High Value) AX: High priority, min safety stock AY: Monitor demand, moderate safety stock AZ: Risky! High carrying cost, high stockout cost. Consider drop-ship or close-to-demand manufacturing
B (Mid Value) BX: Standard EOQ BY: Seasonal forecasting BZ: Caution, may phase out
C (Low Value) CX: Order frequently in small lots CY: Bulk ordering acceptable CZ: Overstock, low priority

Management strategy: - AX: Tight inventory control, reorder frequently, minimize safety stock - AY, AZ: Demand forecasting critical, collaborative planning with customers - BX, BY: Standard practices - CX, CY, CZ: Lenient, bulk ordering acceptable, high service levels can be relaxed

54.2.1 Code: ABC-XYZ Segmentation

library(tidyverse)

# Generate synthetic 500 SKU dataset
set.seed(4617)
sku_data <- tibble(
  sku_id = paste0("SKU_", 1:500),
  annual_demand = rnorm(500, mean = 5000, sd = 2000) |> pmax(100),
  unit_cost = rlnorm(500, meanlog = log(500), sdlog = 0.8),
  daily_demand_sd = rnorm(500, mean = 200, sd = 80) |> pmax(10)
) |>
  mutate(
    annual_value = annual_demand * unit_cost,
    demand_cv = daily_demand_sd / (annual_demand / 365)
  )

# ABC segmentation by value
total_value <- sum(sku_data$annual_value)
abc_data <- sku_data |>
  arrange(desc(annual_value)) |>
  mutate(
    cumulative_value = cumsum(annual_value),
    cumulative_pct = cumulative_value / total_value,
    abc_class = case_when(
      cumulative_pct <= 0.80 ~ "A",
      cumulative_pct <= 0.95 ~ "B",
      TRUE ~ "C"
    )
  )

# XYZ segmentation by demand CV
abc_xyz <- abc_data |>
  mutate(
    xyz_class = case_when(
      demand_cv < 0.5 ~ "X",
      demand_cv < 1.0 ~ "Y",
      TRUE ~ "Z"
    ),
    segment = paste0(abc_class, xyz_class)
  )

# Summary by segment
segment_summary <- abc_xyz |>
  group_by(segment) |>
  summarise(
    n_skus = n(),
    pct_skus = 100 * n() / nrow(abc_xyz),
    total_value_m = sum(annual_value) / 1e6,
    pct_value = 100 * sum(annual_value) / total_value,
    avg_cv = mean(demand_cv),
    .groups = "drop"
  ) |>
  arrange(segment)

print(segment_summary)

# ABC-XYZ matrix heatmap
library(ggplot2)
matrix_data <- abc_xyz |>
  group_by(abc_class, xyz_class) |>
  summarise(
    sku_count = n(),
    inventory_value_m = sum(annual_value) / 1e6,
    .groups = "drop"
  )

ggplot(matrix_data, aes(x = xyz_class, y = abc_class, fill = sku_count)) +
  geom_tile(colour = "white", size = 1) +
  geom_text(aes(label = paste0(sku_count, "\n(₦", round(inventory_value_m, 0), "m)")),
            colour = "white", size = 3) +
  scale_fill_gradient(low = "lightblue", high = "darkblue") +
  labs(title = "ABC-XYZ Segmentation Matrix",
       x = "Demand Variability (XYZ)", y = "Inventory Value (ABC)",
       fill = "SKU Count") +
  theme_minimal()

# Management summary
cat("\n=== ABC-XYZ SEGMENTATION SUMMARY ===\n")
cat("Total SKUs:", nrow(abc_xyz), "\n")
cat("Total Inventory Value: ₦", format(total_value/1e6, digits=0), "m\n\n")

ax_items <- abc_xyz |> filter(segment %in% c("AX", "AY", "AZ"))
cat("A-Items (Focus area):", nrow(ax_items), "SKUs, ₦",
    format(sum(ax_items$annual_value)/1e6, digits=0), "m\n")
cat("  - AX (Priority):", nrow(filter(ax_items, segment=="AX")), "\n")
cat("  - AY (Monitor):", nrow(filter(ax_items, segment=="AY")), "\n")
cat("  - AZ (Risky):", nrow(filter(ax_items, segment=="AZ")), "\n\n")

c_items <- abc_xyz |> filter(abc_class == "C")
cat("C-Items (Low priority):", nrow(c_items), "SKUs, ₦",
    format(sum(c_items$annual_value)/1e6, digits=0), "m\n")
cat("  - Can relax service levels; bulk ordering acceptable\n")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(4617)

# Generate synthetic 500 SKU dataset
sku_data = pd.DataFrame({
    'sku_id': [f'SKU_{i}' for i in range(1, 501)],
    'annual_demand': np.random.normal(5000, 2000, 500).clip(100),
    'unit_cost': np.random.lognormal(np.log(500), 0.8, 500),
    'daily_demand_sd': np.random.normal(200, 80, 500).clip(10)
})

sku_data['annual_value'] = sku_data['annual_demand'] * sku_data['unit_cost']
sku_data['demand_cv'] = sku_data['daily_demand_sd'] / (sku_data['annual_demand'] / 365)

# ABC segmentation by value
total_value = sku_data['annual_value'].sum()
sku_data = sku_data.sort_values('annual_value', ascending=False).reset_index(drop=True)
sku_data['cumulative_value'] = sku_data['annual_value'].cumsum()
sku_data['cumulative_pct'] = sku_data['cumulative_value'] / total_value

sku_data['abc_class'] = pd.cut(
    sku_data['cumulative_pct'],
    bins=[0, 0.80, 0.95, 1.0],
    labels=['A', 'B', 'C'],
    include_lowest=True
)

# XYZ segmentation by demand CV
sku_data['xyz_class'] = pd.cut(
    sku_data['demand_cv'],
    bins=[0, 0.5, 1.0, np.inf],
    labels=['X', 'Y', 'Z'],
    include_lowest=True
)

sku_data['segment'] = sku_data['abc_class'].astype(str) + sku_data['xyz_class'].astype(str)

# Segment summary
segment_summary = sku_data.groupby('segment').agg({
    'sku_id': 'count',
    'annual_value': 'sum',
    'demand_cv': 'mean'
}).rename(columns={'sku_id': 'n_skus', 'annual_value': 'total_value_m'})

segment_summary['pct_skus'] = 100 * segment_summary['n_skus'] / len(sku_data)
segment_summary['pct_value'] = 100 * segment_summary['total_value_m'] / total_value
segment_summary['total_value_m'] = segment_summary['total_value_m'] / 1e6

print(segment_summary.round(2))

# ABC-XYZ matrix heatmap
matrix_data = sku_data.groupby(['abc_class', 'xyz_class']).agg({
    'sku_id': 'count',
    'annual_value': 'sum'
}).rename(columns={'sku_id': 'sku_count', 'annual_value': 'inventory_value'})
matrix_data['inventory_value_m'] = matrix_data['inventory_value'] / 1e6

pivot_data = matrix_data['sku_count'].unstack(fill_value=0)

plt.figure(figsize=(10, 6))
sns.heatmap(pivot_data, annot=True, fmt='d', cmap='Blues', cbar_kws={'label': 'SKU Count'})
plt.title('ABC-XYZ Segmentation Matrix')
plt.xlabel('Demand Variability (XYZ)')
plt.ylabel('Inventory Value (ABC)')
plt.tight_layout()
plt.show()

# Management summary
print(f"\n=== ABC-XYZ SEGMENTATION SUMMARY ===")
print(f"Total SKUs: {len(sku_data)}")
print(f"Total Inventory Value: ₦{total_value/1e6:.0f}m\n")

ax_items = sku_data[sku_data['segment'].isin(['AX', 'AY', 'AZ'])]
print(f"A-Items (Focus area): {len(ax_items)} SKUs, ₦{ax_items['annual_value'].sum()/1e6:.0f}m")
print(f"  - AX (Priority): {len(sku_data[sku_data['segment']=='AX'])}")
print(f"  - AY (Monitor): {len(sku_data[sku_data['segment']=='AY'])}")
print(f"  - AZ (Risky): {len(sku_data[sku_data['segment']=='AZ'])}\n")

c_items = sku_data[sku_data['abc_class'] == 'C']
print(f"C-Items (Low priority): {len(c_items)} SKUs, ₦{c_items['annual_value'].sum()/1e6:.0f}m")
Caution📝 Section 49.2 Review Questions

  1. A company has 500 SKUs with total annual inventory value of ₦500m. ABC analysis shows 50 A-items worth ₦400m, 150 B-items worth ₦80m, 300 C-items worth ₦20m. What percentage of SKUs are A-items? What percentage of value do they represent?

  2. XYZ segmentation measures demand variability via coefficient of variation (CV). Define CV. Why is CV > 1 considered erratic, and how should inventory managers treat Z items differently from X items?

  3. A product has mean daily demand = 100 units with SD = 80 units. Calculate its CV. Is it X, Y, or Z? How should the company treat this item?

  4. Explain the rationale for lower service levels (e.g., 80%) for C items versus higher service levels (e.g., 98%) for A items. What tradeoff is being made?

  5. An AZ item is high-value but erratic in demand. Why is this “risky”? What three alternative strategies could a manager use to reduce risk?

54.3 Economic Order Quantity (EOQ) and Sensitivity

The EOQ formula \(\sqrt{\frac{2Ds}{h}}\) balances holding and ordering costs. For a pharmaceutical SKU:

  • D = 10,000 units/year
  • s = ₦100,000 per order (administrative + shipping)
  • h = ₦50 per unit-year (carrying cost)

\[\text{EOQ} = \sqrt{\frac{2 \times 10,000 \times 100,000}{50}} = \sqrt{40,000,000} = 6,325 \text{ units}\]

Order 6,325 units roughly every 231 days (10,000 / 6,325 × 365). Average inventory = EOQ/2 = 3,162 units. Annual holding cost = 3,162 × 50 = ₦158k. Annual ordering cost = (10,000 / 6,325) × 100k = ₦158k. Total = ₦316k.

Sensitivity analysis examines how EOQ changes with parameter shifts:

  • If holding cost doubles (₦100/unit), EOQ drops to 4,472 units (27% decrease) → order more frequently, carry less inventory
  • If ordering cost halves (₦50k), EOQ drops to 4,472 units → frequent small orders become economical
  • If demand doubles (20,000 units/year), EOQ rises to 8,944 units (41% increase)

EOQ is insensitive to small parameter changes; a ±10% swing in parameters changes EOQ by ~5%. This robustness makes EOQ practical despite real-world uncertainty.

Note📘 Theory: EOQ Assumptions and Extensions

EOQ assumes: (1) constant demand, (2) instantaneous replenishment, (3) no quantity discounts, (4) no stockouts. Real supply chains violate all assumptions. Extensions: variable demand (Stochastic EOQ), lead time (reorder point model), quantity discounts (tiered EOQ), and multiple items with constraints (linear programming). For most A items, the stochastic EOQ (accounting for demand variability) is more realistic.

54.3.1 Code: EOQ Calculation and Sensitivity Analysis

# EOQ function
eoq <- function(D, s, h) {
  sqrt(2 * D * s / h)
}

# Total cost function
total_inventory_cost <- function(Q, D, s, h) {
  (Q / 2) * h + (D / Q) * s
}

# Example: Pharmaceutical SKU
D <- 10000  # annual demand (units)
s <- 100000 # ordering cost per order (₦)
h <- 50     # holding cost per unit-year (₦)

# Calculate EOQ
eoq_qty <- eoq(D, s, h)
orders_per_year <- D / eoq_qty
avg_inventory <- eoq_qty / 2
annual_holding_cost <- avg_inventory * h
annual_ordering_cost <- orders_per_year * s
total_cost <- annual_holding_cost + annual_ordering_cost

cat("=== EOQ CALCULATION ===\n")
cat(sprintf("Annual Demand (D): %.0f units\n", D))
cat(sprintf("Ordering Cost (s): ₦%.0f per order\n", s))
cat(sprintf("Holding Cost (h): ₦%.2f per unit-year\n", h))
cat(sprintf("\nEOQ: %.0f units\n", eoq_qty))
cat(sprintf("Orders/Year: %.2f\n", orders_per_year))
cat(sprintf("Days between orders: %.0f days\n", 365 / orders_per_year))
cat(sprintf("Avg Inventory: %.0f units\n", avg_inventory))
cat(sprintf("Annual Holding Cost: ₦%.0f\n", annual_holding_cost))
cat(sprintf("Annual Ordering Cost: ₦%.0f\n", annual_ordering_cost))
cat(sprintf("Total Annual Cost: ₦%.0f\n\n", total_cost))

# Sensitivity analysis: vary holding cost and ordering cost ±50%
library(tidyverse)

h_range <- h * c(0.5, 0.75, 1, 1.25, 1.5)
s_range <- s * c(0.5, 0.75, 1, 1.25, 1.5)

sensitivity <- expand_grid(h_var = h_range, s_var = s_range) |>
  mutate(
    eoq_var = eoq(D, s_var, h_var),
    pct_change_eoq = 100 * (eoq_var - eoq_qty) / eoq_qty,
    total_cost_var = total_inventory_cost(eoq_var, D, s_var, h_var)
  )

cat("=== SENSITIVITY ANALYSIS: EOQ Impact ===\n")
cat("(% change in EOQ vs base case)\n\n")

# Holding cost sensitivity
h_sens <- sensitivity |>
  filter(s_var == s) |>
  select(h_var, eoq_var, pct_change_eoq) |>
  mutate(h_pct_change = 100 * (h_var - h) / h) |>
  arrange(h_var)

cat("Holding Cost Sensitivity:\n")
print(h_sens |> select(h_pct_change, eoq_var, pct_change_eoq) |>
  rename("h % change" = h_pct_change, "EOQ" = eoq_var, "EOQ % change" = pct_change_eoq))

cat("\nOrdering Cost Sensitivity:\n")

s_sens <- sensitivity |>
  filter(h_var == h) |>
  select(s_var, eoq_var, pct_change_eoq) |>
  mutate(s_pct_change = 100 * (s_var - s) / s) |>
  arrange(s_var)

print(s_sens |> select(s_pct_change, eoq_var, pct_change_eoq) |>
  rename("s % change" = s_pct_change, "EOQ" = eoq_var, "EOQ % change" = pct_change_eoq))

# Heatmap: EOQ values across holding and ordering cost ranges
library(ggplot2)

heatmap_data <- sensitivity |>
  mutate(
    h_label = paste0(round(100*(h_var-h)/h, 0), "%"),
    s_label = paste0(round(100*(s_var-s)/s, 0), "%")
  )

ggplot(heatmap_data, aes(x = s_label, y = h_label, fill = eoq_var)) +
  geom_tile(colour = "white") +
  geom_text(aes(label = round(eoq_var)), size = 2.5) +
  scale_fill_gradient(low = "lightgreen", high = "darkgreen") +
  labs(title = "EOQ Sensitivity: Holding Cost vs Ordering Cost",
       x = "Ordering Cost Change (%)", y = "Holding Cost Change (%)",
       fill = "EOQ (units)") +
  theme_minimal()

# Cost comparison at different order quantities
quantities <- seq(2000, 12000, by = 500)
cost_comparison <- tibble(
  quantity = quantities,
  cost_at_qty = sapply(quantities, function(q) total_inventory_cost(q, D, s, h)),
  cost_at_eoq = total_cost
)

ggplot(cost_comparison, aes(x = quantity, y = cost_at_qty)) +
  geom_line(linewidth = 1) +
  geom_point(data = tibble(quantity = eoq_qty, cost_at_qty = total_cost),
             colour = "red", size = 3) +
  geom_vline(xintercept = eoq_qty, linetype = "dashed", colour = "red") +
  annotate("text", x = eoq_qty, y = total_cost * 1.05,
           label = paste0("EOQ = ", format(round(eoq_qty), big.mark=",")),
           hjust = 0.5) +
  labs(title = "Total Inventory Cost vs Order Quantity",
       x = "Order Quantity (units)", y = "Annual Cost (₦)") +
  theme_minimal()
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize_scalar

def eoq(D, s, h):
    """Calculate Economic Order Quantity"""
    return np.sqrt(2 * D * s / h)

def total_inventory_cost(Q, D, s, h):
    """Calculate total inventory cost for a given order quantity"""
    return (Q / 2) * h + (D / Q) * s

# Example: Pharmaceutical SKU
D = 10000   # annual demand (units)
s = 100000  # ordering cost per order (₦)
h = 50      # holding cost per unit-year (₦)

# Calculate EOQ
eoq_qty = eoq(D, s, h)
orders_per_year = D / eoq_qty
avg_inventory = eoq_qty / 2
annual_holding_cost = avg_inventory * h
annual_ordering_cost = orders_per_year * s
total_cost = annual_holding_cost + annual_ordering_cost

print("=== EOQ CALCULATION ===")
print(f"Annual Demand (D): {D:,.0f} units")
print(f"Ordering Cost (s): ₦{s:,.0f} per order")
print(f"Holding Cost (h): ₦{h:.2f} per unit-year")
print(f"\nEOQ: {eoq_qty:,.0f} units")
print(f"Orders/Year: {orders_per_year:.2f}")
print(f"Days between orders: {365/orders_per_year:.0f} days")
print(f"Avg Inventory: {avg_inventory:,.0f} units")
print(f"Annual Holding Cost: ₦{annual_holding_cost:,.0f}")
print(f"Annual Ordering Cost: ₦{annual_ordering_cost:,.0f}")
print(f"Total Annual Cost: ₦{total_cost:,.0f}\n")

# Sensitivity analysis: vary holding cost and ordering cost ±50%
h_range = h * np.array([0.5, 0.75, 1, 1.25, 1.5])
s_range = s * np.array([0.5, 0.75, 1, 1.25, 1.5])

# Holding cost sensitivity
h_sensitivity = []
for h_var in h_range:
    eoq_var = eoq(D, s, h_var)
    pct_change = 100 * (eoq_var - eoq_qty) / eoq_qty
    h_pct_change = 100 * (h_var - h) / h
    h_sensitivity.append({
        'h_pct_change': h_pct_change,
        'eoq': eoq_var,
        'eoq_pct_change': pct_change
    })

h_sens_df = pd.DataFrame(h_sensitivity)
print("=== SENSITIVITY ANALYSIS: EOQ Impact ===")
print("Holding Cost Sensitivity:\n")
print(h_sens_df.round(2))

# Ordering cost sensitivity
s_sensitivity = []
for s_var in s_range:
    eoq_var = eoq(D, s_var, h)
    pct_change = 100 * (eoq_var - eoq_qty) / eoq_qty
    s_pct_change = 100 * (s_var - s) / s
    s_sensitivity.append({
        's_pct_change': s_pct_change,
        'eoq': eoq_var,
        'eoq_pct_change': pct_change
    })

s_sens_df = pd.DataFrame(s_sensitivity)
print("\nOrdering Cost Sensitivity:\n")
print(s_sens_df.round(2))

# Heatmap: EOQ values across parameter ranges
sensitivity_matrix = np.zeros((len(h_range), len(s_range)))
for i, h_var in enumerate(h_range):
    for j, s_var in enumerate(s_range):
        sensitivity_matrix[i, j] = eoq(D, s_var, h_var)

fig, ax = plt.subplots(figsize=(10, 6))
sns.heatmap(sensitivity_matrix, annot=True, fmt='.0f', cmap='Greens',
            xticklabels=[f'{100*(x-s)/s:+.0f}%' for x in s_range],
            yticklabels=[f'{100*(x-h)/h:+.0f}%' for x in h_range],
            cbar_kws={'label': 'EOQ (units)'})
ax.set_xlabel('Ordering Cost Change (%)')
ax.set_ylabel('Holding Cost Change (%)')
ax.set_title('EOQ Sensitivity: Holding Cost vs Ordering Cost')
plt.tight_layout()
plt.show()

# Cost comparison at different order quantities
quantities = np.arange(2000, 12001, 500)
costs = [total_inventory_cost(q, D, s, h) for q in quantities]

plt.figure(figsize=(10, 6))
plt.plot(quantities, costs, linewidth=2, label='Total Cost')
plt.plot(eoq_qty, total_cost, 'ro', markersize=8, label=f'EOQ = {eoq_qty:,.0f}')
plt.axvline(eoq_qty, linestyle='--', color='red', alpha=0.5)
plt.xlabel('Order Quantity (units)')
plt.ylabel('Annual Cost (₦)')
plt.title('Total Inventory Cost vs Order Quantity')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Caution📝 Section 49.3 Review Questions

  1. State the EOQ formula and its four key assumptions. Which assumptions are most frequently violated in practice?

  2. Annual demand = 12,000 units, ordering cost = ₦80,000, holding cost = ₦40/unit-year. Calculate EOQ, number of orders per year, average inventory, and total annual cost.

  3. If holding cost increases from ₦40 to ₦60 per unit-year (50% increase), how does EOQ change? By what percentage?

  4. Explain why EOQ is said to be “robust” to parameter uncertainty. If your estimate of ordering cost is off by 20%, what is the maximum EOQ error?

  5. A manager argues: “We should order large quantities infrequently because ordering cost is high.” Under what condition (in terms of holding cost) would this be optimal?

54.4 Reorder Point and Safety Stock

In reality, demand and lead time vary. A simple reorder point (RP) approach:

\[\text{RP} = \text{Average Daily Demand} \times \text{Average Lead Time} + \text{Safety Stock}\]

Safety stock buffers against both demand spikes and delayed replenishment. If daily demand averages 30 units, lead time 10 days, and 95% service level required:

Average demand during lead time = 30 × 10 = 300 units. But demand is uncertain; a z-score of 1.645 covers 95% of the normal distribution. If daily demand standard deviation = 5 units, then lead-time demand SD = 5 × √10 ≈ 15.8 units.

\[\text{Safety Stock} = z \times \sigma_{\text{LT}} = 1.645 \times 15.8 ≈ 26 \text{ units}\]

\[\text{RP} = 300 + 26 = 326 \text{ units}\]

Reorder when inventory hits 326 units. Combining RP with EOQ: order 6,325 units when inventory reaches 326 units.

Higher service levels require higher safety stock (and higher carrying cost). Service level = 90% → z = 1.28, safety stock ≈ 20 units. Service level = 99% → z = 2.33, safety stock ≈ 37 units. The company chooses service level based on stockout cost vs holding cost.

Note📘 Theory: Service Level and Stockout Probability

Service level = P(Demand ≤ Inventory available during lead time). For normal demand, service level corresponds to a z-score: 95% service = z=1.645, 99% service = z=2.33. The optimal service level equates marginal cost of safety stock (holding cost per unit) to marginal benefit (avoided stockout cost per unit shortage).

54.4.1 Code: Safety Stock Calculation

library(tidyverse)

# Service level lookup table
service_level_z <- tibble(
  service_level_pct = c(90, 95, 99),
  z_score = c(1.282, 1.645, 2.326)
)

# Safety stock function
safety_stock <- function(z, sigma_d, mean_d, sigma_L, mean_L) {
  # Demand during lead time variance: Var = sigma_d^2 * L + mean_d^2 * sigma_L^2
  var_lt_demand <- (sigma_d^2 * mean_L) + (mean_d^2 * sigma_L^2)
  sigma_lt <- sqrt(var_lt_demand)
  ss <- z * sigma_lt
  return(ss)
}

# Reorder point function
reorder_point <- function(mean_d, mean_L, ss) {
  return(mean_d * mean_L + ss)
}

# Example 1: High-demand stable product
cat("=== EXAMPLE 1: High-Demand Stable Product ===\n")
mean_d_1 <- 100      # daily demand mean (units)
sigma_d_1 <- 12      # daily demand SD (units)
mean_L_1 <- 12       # lead time mean (days)
sigma_L_1 <- 1.5     # lead time SD (days)

# Calculate for 95% service level
z_95 <- 1.645
ss_95 <- safety_stock(z_95, sigma_d_1, mean_d_1, sigma_L_1, mean_L_1)
rp_95 <- reorder_point(mean_d_1, mean_L_1, ss_95)

cat(sprintf("Daily demand: μ=%.0f, σ=%.0f units\n", mean_d_1, sigma_d_1))
cat(sprintf("Lead time: μ=%.0f, σ=%.1f days\n\n", mean_L_1, sigma_L_1))

cat("95% Service Level:\n")
cat(sprintf("  Safety Stock: %.0f units\n", ss_95))
cat(sprintf("  Reorder Point: %.0f units\n\n", rp_95))

# Calculate for 99% service level
z_99 <- 2.326
ss_99 <- safety_stock(z_99, sigma_d_1, mean_d_1, sigma_L_1, mean_L_1)
rp_99 <- reorder_point(mean_d_1, mean_L_1, ss_99)

cat("99% Service Level:\n")
cat(sprintf("  Safety Stock: %.0f units\n", ss_99))
cat(sprintf("  Reorder Point: %.0f units\n", rp_99))

# Cost comparison
h <- 50              # holding cost per unit-year (₦)
annual_ss_cost_95 <- ss_95 * h
annual_ss_cost_99 <- ss_99 * h
cost_increase <- annual_ss_cost_99 - annual_ss_cost_95

cat(sprintf("\nAnnual Holding Cost (Safety Stock):\n"))
cat(sprintf("  95% SL: ₦%.0f\n", annual_ss_cost_95))
cat(sprintf("  99% SL: ₦%.0f\n", annual_ss_cost_99))
cat(sprintf("  Increase: ₦%.0f (%.1f%% higher)\n\n", cost_increase,
    100*cost_increase/annual_ss_cost_95))

# Example 2: Erratic demand product
cat("=== EXAMPLE 2: Erratic Demand Product ===\n")
mean_d_2 <- 25       # daily demand mean (units)
sigma_d_2 <- 22      # daily demand SD (units) - high variability
mean_L_2 <- 18       # lead time mean (days)
sigma_L_2 <- 2.5     # lead time SD (days)

ss_95_2 <- safety_stock(z_95, sigma_d_2, mean_d_2, sigma_L_2, mean_L_2)
rp_95_2 <- reorder_point(mean_d_2, mean_L_2, ss_95_2)
ss_99_2 <- safety_stock(z_99, sigma_d_2, mean_d_2, sigma_L_2, mean_L_2)
rp_99_2 <- reorder_point(mean_d_2, mean_L_2, ss_99_2)

cat(sprintf("Daily demand: μ=%.0f, σ=%.0f units (CV=%.2f - ERRATIC)\n",
    mean_d_2, sigma_d_2, sigma_d_2/mean_d_2))
cat(sprintf("Lead time: μ=%.0f, σ=%.1f days\n\n", mean_L_2, sigma_L_2))

cat("95% Service Level:\n")
cat(sprintf("  Safety Stock: %.0f units\n", ss_95_2))
cat(sprintf("  Reorder Point: %.0f units\n\n", rp_95_2))

cat("99% Service Level:\n")
cat(sprintf("  Safety Stock: %.0f units\n", ss_99_2))
cat(sprintf("  Reorder Point: %.0f units\n", rp_99_2))

# Service level vs safety stock curve
service_levels <- seq(80, 99, by = 1)
z_scores <- qnorm(service_levels / 100)
safety_stocks <- sapply(z_scores, function(z) safety_stock(z, sigma_d_1, mean_d_1, sigma_L_1, mean_L_1))
ss_costs <- safety_stocks * h

library(ggplot2)

ggplot(tibble(sl = service_levels, ss = safety_stocks, cost = ss_costs),
       aes(x = sl, y = ss)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(title = "Safety Stock vs Service Level",
       x = "Service Level (%)", y = "Safety Stock (units)") +
  theme_minimal()

ggplot(tibble(sl = service_levels, cost = ss_costs),
       aes(x = sl, y = cost)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(title = "Annual Holding Cost (Safety Stock) vs Service Level",
       x = "Service Level (%)", y = "Annual Cost (₦)") +
  theme_minimal()
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm

def safety_stock(z, sigma_d, mean_d, sigma_L, mean_L):
    """Calculate safety stock for given service level"""
    var_lt_demand = (sigma_d**2 * mean_L) + (mean_d**2 * sigma_L**2)
    sigma_lt = np.sqrt(var_lt_demand)
    ss = z * sigma_lt
    return ss

def reorder_point(mean_d, mean_L, ss):
    """Calculate reorder point"""
    return mean_d * mean_L + ss

# Service level lookup
service_level_z = {
    90: 1.282,
    95: 1.645,
    99: 2.326
}

# Example 1: High-demand stable product
print("=== EXAMPLE 1: High-Demand Stable Product ===")
mean_d_1 = 100       # daily demand mean (units)
sigma_d_1 = 12       # daily demand SD (units)
mean_L_1 = 12        # lead time mean (days)
sigma_L_1 = 1.5      # lead time SD (days)

z_95 = 1.645
ss_95 = safety_stock(z_95, sigma_d_1, mean_d_1, sigma_L_1, mean_L_1)
rp_95 = reorder_point(mean_d_1, mean_L_1, ss_95)

z_99 = 2.326
ss_99 = safety_stock(z_99, sigma_d_1, mean_d_1, sigma_L_1, mean_L_1)
rp_99 = reorder_point(mean_d_1, mean_L_1, ss_99)

print(f"Daily demand: μ={mean_d_1:.0f}, σ={sigma_d_1:.0f} units")
print(f"Lead time: μ={mean_L_1:.0f}, σ={sigma_L_1:.1f} days\n")
print(f"95% Service Level:")
print(f"  Safety Stock: {ss_95:.0f} units")
print(f"  Reorder Point: {rp_95:.0f} units\n")
print(f"99% Service Level:")
print(f"  Safety Stock: {ss_99:.0f} units")
print(f"  Reorder Point: {rp_99:.0f} units\n")

# Cost comparison
h = 50  # holding cost per unit-year (₦)
annual_ss_cost_95 = ss_95 * h
annual_ss_cost_99 = ss_99 * h
cost_increase = annual_ss_cost_99 - annual_ss_cost_95

print(f"Annual Holding Cost (Safety Stock):")
print(f"  95% SL: ₦{annual_ss_cost_95:.0f}")
print(f"  99% SL: ₦{annual_ss_cost_99:.0f}")
print(f"  Increase: ₦{cost_increase:.0f} ({100*cost_increase/annual_ss_cost_95:.1f}% higher)\n")

# Example 2: Erratic demand product
print("=== EXAMPLE 2: Erratic Demand Product ===")
mean_d_2 = 25        # daily demand mean (units)
sigma_d_2 = 22       # daily demand SD (units)
mean_L_2 = 18        # lead time mean (days)
sigma_L_2 = 2.5      # lead time SD (days)

ss_95_2 = safety_stock(z_95, sigma_d_2, mean_d_2, sigma_L_2, mean_L_2)
rp_95_2 = reorder_point(mean_d_2, mean_L_2, ss_95_2)
ss_99_2 = safety_stock(z_99, sigma_d_2, mean_d_2, sigma_L_2, mean_L_2)
rp_99_2 = reorder_point(mean_d_2, mean_L_2, ss_99_2)

print(f"Daily demand: μ={mean_d_2:.0f}, σ={sigma_d_2:.0f} units (CV={sigma_d_2/mean_d_2:.2f} - ERRATIC)")
print(f"Lead time: μ={mean_L_2:.0f}, σ={sigma_L_2:.1f} days\n")
print(f"95% Service Level:")
print(f"  Safety Stock: {ss_95_2:.0f} units")
print(f"  Reorder Point: {rp_95_2:.0f} units\n")
print(f"99% Service Level:")
print(f"  Safety Stock: {ss_99_2:.0f} units")
print(f"  Reorder Point: {rp_99_2:.0f} units\n")

# Service level vs safety stock curve
service_levels = np.arange(80, 100, 1)
z_scores = norm.ppf(service_levels / 100)
safety_stocks = [safety_stock(z, sigma_d_1, mean_d_1, sigma_L_1, mean_L_1) for z in z_scores]
ss_costs = [ss * h for ss in safety_stocks]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

ax1.plot(service_levels, safety_stocks, linewidth=2, marker='o')
ax1.set_xlabel('Service Level (%)')
ax1.set_ylabel('Safety Stock (units)')
ax1.set_title('Safety Stock vs Service Level')
ax1.grid(True, alpha=0.3)

ax2.plot(service_levels, ss_costs, linewidth=2, marker='o', color='orange')
ax2.set_xlabel('Service Level (%)')
ax2.set_ylabel('Annual Cost (₦)')
ax2.set_title('Annual Holding Cost (Safety Stock) vs Service Level')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
Caution📝 Section 49.4 Review Questions

  1. Define reorder point (RP) and safety stock (SS). Write the formula for RP in terms of demand, lead time, and SS.

  2. Daily demand is normally distributed with mean = 50 units and SD = 8 units. Lead time is constant at 14 days. Calculate safety stock for 95% and 99% service levels. What is the reorder point for each?

  3. A product has variable demand AND variable lead time. The formula for lead-time demand variance is \(\text{Var} = \sigma_d^2 \times L + \mu_d^2 \times \sigma_L^2\). Explain each term: what does the first term represent? The second?

  4. Increasing service level from 95% to 99% increases safety stock cost. Under what condition (in terms of stockout cost) is 99% justified over 95%?

  5. An ABC analysis shows that an AX item generates ₦2m margin per year but has erratic demand (Z class). Should this item have a higher or lower service level than a BX item? Explain your reasoning.

54.5 Linear Programming for Supply Chain Optimisation

The transportation problem optimises shipping from M warehouses to N customers to minimize total cost subject to supply and demand constraints.

Minimise: \(\sum_i \sum_j c_{ij} x_{ij}\) (total shipping cost)

Subject to: - \(\sum_j x_{ij} = S_i\) (warehouse \(i\) supplies \(S_i\) units) - \(\sum_i x_{ij} = D_j\) (customer \(j\) demands \(D_j\) units) - \(x_{ij} \geq 0\)

where \(c_{ij}\) = cost to ship one unit from warehouse \(i\) to customer \(j\), and \(x_{ij}\) = quantity shipped.

This is a linear program solvable in polynomial time (simplex or network algorithms). For a pan-African distributor with 3 warehouses (Lagos, Abuja, Kano) and 6 customer zones, LP finds the minimum-cost allocation.

Solutions provide: (1) optimal allocation (which warehouse supplies which customer), (2) shadow prices (cost of relaxing a constraint), (3) sensitivity ranges (within what parameter ranges is the solution optimal).

54.5.1 Code: LP Transportation Problem

library(lpSolve)

# Transportation problem setup
# 3 warehouses: Lagos, Abuja, Kano
# 4 customer zones: SW, NW, NC, NE

# Supply from each warehouse (in units, e.g., millions of naira worth)
supply <- c(250, 150, 100)  # Lagos, Abuja, Kano
warehouse_names <- c("Lagos", "Abuja", "Kano")

# Demand from each zone
demand <- c(200, 120, 100, 80)  # SW, NW, NC, NE
zone_names <- c("SW", "NW", "NC", "NE")

# Shipping cost matrix (3x4): cost per unit from warehouse i to zone j
# In ₦ per unit (e.g., ₦100 per unit from Lagos to SW)
cost_matrix <- matrix(c(
  100, 150, 180, 220,   # Lagos to SW, NW, NC, NE
  160, 110, 140, 180,   # Abuja to SW, NW, NC, NE
  200, 170, 120, 80     # Kano to SW, NW, NC, NE
), nrow = 3, byrow = TRUE)

colnames(cost_matrix) <- zone_names
rownames(cost_matrix) <- warehouse_names

cat("=== TRANSPORTATION PROBLEM ===\n")
cat("Supply from warehouses:\n")
print(tibble(Warehouse = warehouse_names, Supply = supply))
cat("\nDemand from zones:\n")
print(tibble(Zone = zone_names, Demand = demand))
cat("\nShipping cost matrix (₦ per unit):\n")
print(cost_matrix)

# Flatten cost matrix for lp.transport
c_vector <- as.vector(t(cost_matrix))

# Solve transportation problem using lpSolve
result <- lp.transport(
  cost.mat = cost_matrix,
  row.signs = rep("=", 3),
  row.rhs = supply,
  col.signs = rep("=", 4),
  col.rhs = demand
)

cat("\n=== OPTIMAL SOLUTION ===\n")
cat(sprintf("Minimum Total Cost: ₦%.0f\n\n", result$objval))

# Extract optimal allocation
optimal_allocation <- matrix(result$solution, nrow = 3, byrow = TRUE)
colnames(optimal_allocation) <- zone_names
rownames(optimal_allocation) <- warehouse_names

cat("Optimal Allocation (units):\n")
print(as.data.frame(optimal_allocation))

cat("\nOptimal Allocation Summary:\n")
allocation_df <- expand_grid(
  Warehouse = warehouse_names,
  Zone = zone_names
) |>
  mutate(
    Units = as.vector(t(optimal_allocation)),
    Cost_per_Unit = as.vector(t(cost_matrix)),
    Total_Cost = Units * Cost_per_Unit
  ) |>
  filter(Units > 0)

print(allocation_df)

cat(sprintf("\nTotal Cost Breakdown:\n"))
print(allocation_df |>
  group_by(Warehouse) |>
  summarise(Total_Cost = sum(Total_Cost), .groups = "drop") |>
  arrange(desc(Total_Cost)))

# Sensitivity analysis: what if shipping cost changes?
cost_increase_pct <- 10  # 10% increase
cost_matrix_sensitive <- cost_matrix * (1 + cost_increase_pct / 100)

result_sensitive <- lp.transport(
  cost.mat = cost_matrix_sensitive,
  row.signs = rep("=", 3),
  row.rhs = supply,
  col.signs = rep("=", 4),
  col.rhs = demand
)

cat(sprintf("\nSensitivity: +%d%% Shipping Cost\n", cost_increase_pct))
cat(sprintf("New Minimum Total Cost: ₦%.0f\n", result_sensitive$objval))
cat(sprintf("Cost Increase: ₦%.0f (%.2f%%)\n",
    result_sensitive$objval - result$objval,
    100 * (result_sensitive$objval - result$objval) / result$objval))
import numpy as np
import pandas as pd
from pulp import *

# Transportation problem setup
warehouse_names = ["Lagos", "Abuja", "Kano"]
zone_names = ["SW", "NW", "NC", "NE"]

# Supply and demand
supply = [250, 150, 100]
demand = [200, 120, 100, 80]

# Shipping cost matrix
cost_matrix = np.array([
    [100, 150, 180, 220],   # Lagos
    [160, 110, 140, 180],   # Abuja
    [200, 170, 120, 80]     # Kano
])

print("=== TRANSPORTATION PROBLEM ===")
print("\nSupply from warehouses:")
print(pd.DataFrame({'Warehouse': warehouse_names, 'Supply': supply}))

print("\nDemand from zones:")
print(pd.DataFrame({'Zone': zone_names, 'Demand': demand}))

print("\nShipping cost matrix (₦ per unit):")
print(pd.DataFrame(cost_matrix, index=warehouse_names, columns=zone_names))

# Create the LP problem
prob = LpProblem("Transportation_Problem", LpMinimize)

# Decision variables: x[i][j] = units shipped from warehouse i to zone j
x = {}
for i in range(len(warehouse_names)):
    for j in range(len(zone_names)):
        x[(i, j)] = LpVariable(f"x_{warehouse_names[i]}_{zone_names[j]}", lowBound=0)

# Objective function: minimize total shipping cost
prob += lpSum([cost_matrix[i, j] * x[(i, j)] for i in range(3) for j in range(4)])

# Supply constraints: each warehouse supplies exactly its supply amount
for i in range(len(warehouse_names)):
    prob += lpSum([x[(i, j)] for j in range(4)]) == supply[i], f"Supply_{warehouse_names[i]}"

# Demand constraints: each zone receives exactly its demand
for j in range(len(zone_names)):
    prob += lpSum([x[(i, j)] for i in range(3)]) == demand[j], f"Demand_{zone_names[j]}"

# Solve the problem
prob.solve(PULP_CBC_CMD(msg=0))

print("\n=== OPTIMAL SOLUTION ===")
print(f"Minimum Total Cost: ₦{value(prob.objective):,.0f}")

# Extract optimal allocation
optimal_allocation = np.zeros((3, 4))
for i in range(3):
    for j in range(4):
        optimal_allocation[i, j] = value(x[(i, j)])

print("\nOptimal Allocation (units):")
print(pd.DataFrame(optimal_allocation, index=warehouse_names, columns=zone_names))

# Allocation summary
allocation_list = []
for i in range(3):
    for j in range(4):
        units = value(x[(i, j)])
        if units > 0:
            allocation_list.append({
                'Warehouse': warehouse_names[i],
                'Zone': zone_names[j],
                'Units': units,
                'Cost_per_Unit': cost_matrix[i, j],
                'Total_Cost': units * cost_matrix[i, j]
            })

allocation_df = pd.DataFrame(allocation_list)
print("\nOptimal Allocation Summary:")
print(allocation_df.to_string(index=False))

print("\nTotal Cost by Warehouse:")
print(allocation_df.groupby('Warehouse')['Total_Cost'].sum().sort_values(ascending=False))

# Sensitivity analysis
print("\n=== SENSITIVITY ANALYSIS ===")
cost_increase_pct = 10
cost_matrix_sensitive = cost_matrix * (1 + cost_increase_pct / 100)

prob_sensitive = LpProblem("Transportation_Sensitive", LpMinimize)
x_sensitive = {}
for i in range(3):
    for j in range(4):
        x_sensitive[(i, j)] = LpVariable(f"x_{i}_{j}", lowBound=0)

prob_sensitive += lpSum([cost_matrix_sensitive[i, j] * x_sensitive[(i, j)]
                          for i in range(3) for j in range(4)])

for i in range(3):
    prob_sensitive += lpSum([x_sensitive[(i, j)] for j in range(4)]) == supply[i]

for j in range(4):
    prob_sensitive += lpSum([x_sensitive[(i, j)] for i in range(3)]) == demand[j]

prob_sensitive.solve(PULP_CBC_CMD(msg=0))

new_cost = value(prob_sensitive.objective)
cost_increase_amt = new_cost - value(prob.objective)
print(f"After +{cost_increase_pct}% Shipping Cost Increase:")
print(f"New Minimum Total Cost: ₦{new_cost:,.0f}")
print(f"Cost Increase: ₦{cost_increase_amt:,.0f} ({100*cost_increase_amt/value(prob.objective):.2f}%)")
Caution📝 Section 49.5 Review Questions

  1. In the transportation problem, what do the constraints \(\sum_j x_{ij} = S_i\) and \(\sum_i x_{ij} = D_j\) represent? Why are both necessary?

  2. A warehouse has excess supply (supply > total demand). How would you modify the LP to allow unshipped inventory? What variable would you add?

  3. Explain what “shadow price” means in the context of the transportation problem. For the demand constraint for zone SW, what does a shadow price of ₦50 indicate?

  4. The transportation solution should allocate 3 units from Lagos to NW and 4 units from Abuja to NW, but your code only returns one route. Why might this happen, and is it correct?

  5. How would you extend the transportation problem to include warehouse capacity constraints (e.g., Lagos warehouse can hold max ₦250m worth of inventory but may not use all of it)?

54.6 Simulation for Inventory Management

Monte Carlo simulation models stochastic inventory policies and compares them. Simulate 365 days of demand (draw from historical distribution, add seasonality). For each day:

  1. Demand \(d_t\) arrives
  2. Fulfil from inventory; update stock
  3. If stock < reorder point, place order for Q units
  4. Order arrives after lead time; update inventory
  5. Record stockout (if inventory < demand) and carrying cost

Compare two policies: (s, Q) (order Q units when inventory falls below s) vs (s, S) (order up to S when inventory falls below s). Metrics: average cost (holding + ordering + stockout), service level achieved, maximum inventory level needed.

Simulation is valuable when closed-form solutions are intractable (e.g., nonlinear demand, multiple constraints, complex lead time distributions).

54.6.1 Code: Monte Carlo Inventory Simulation

library(tidyverse)

# Simulation parameters
set.seed(8293)
n_simulations <- 1000
horizon <- 365  # days

# Inventory policy parameters
mean_daily_demand <- 50
sd_daily_demand <- 8
lead_time <- 14  # days
ordering_cost <- 100000  # ₦ per order
holding_cost <- 50  # ₦ per unit-year
stockout_cost <- 500  # ₦ per unit of shortage

# Policy 1: (s, Q) - order Q when inventory drops below s
s_1 <- 500  # reorder point
Q_1 <- 600  # order quantity

# Policy 2: (s, S) - order up to S when inventory drops below s
s_2 <- 500
S_2 <- 1200

# Simulation function for (s, Q) policy
simulate_sq_policy <- function(n_days, s, Q, mean_d, sd_d, L, order_cost, hold_cost, stockout_cost) {
  inventory <- Q  # Start with one order
  on_order <- 0   # Units on order (not yet received)
  days_until_arrival <- 0

  total_holding_cost <- 0
  total_ordering_cost <- 0
  total_stockout_cost <- 0
  num_stockouts <- 0
  max_inventory <- 0
  order_count <- 0

  for (day in 1:n_days) {
    # Check if order arrives
    if (on_order > 0 && days_until_arrival == 0) {
      inventory <- inventory + on_order
      on_order <- 0
    } else if (on_order > 0) {
      days_until_arrival <- days_until_arrival - 1
    }

    # Daily demand
    demand <- max(0, rnorm(1, mean_d, sd_d))

    # Fulfil demand
    if (inventory >= demand) {
      inventory <- inventory - demand
    } else {
      shortage <- demand - inventory
      total_stockout_cost <- total_stockout_cost + shortage * stockout_cost
      num_stockouts <- num_stockouts + 1
      inventory <- 0
    }

    # Holding cost (daily, annualized)
    total_holding_cost <- total_holding_cost + (inventory * hold_cost / 365)

    # Track max inventory
    max_inventory <- max(max_inventory, inventory)

    # Reorder if inventory falls below s
    if (inventory < s && on_order == 0) {
      on_order <- Q
      days_until_arrival <- L
      total_ordering_cost <- total_ordering_cost + order_cost
      order_count <- order_count + 1
    }
  }

  service_level <- 1 - (num_stockouts / n_days)

  return(list(
    total_cost = total_holding_cost + total_ordering_cost + total_stockout_cost,
    holding_cost = total_holding_cost,
    ordering_cost = total_ordering_cost,
    stockout_cost = total_stockout_cost,
    avg_inventory = mean(inventory),  # Note: simple approximation
    max_inventory = max_inventory,
    service_level = service_level,
    num_stockouts = num_stockouts,
    orders = order_count
  ))
}

# Run simulations
cat("=== INVENTORY POLICY SIMULATION ===\n")
cat(sprintf("Simulations: %d runs × %d days\n\n", n_simulations, horizon))

results_sq <- list()
results_ss <- list()

for (sim in 1:n_simulations) {
  results_sq[[sim]] <- simulate_sq_policy(horizon, s_1, Q_1, mean_daily_demand, sd_daily_demand,
                                          lead_time, ordering_cost, holding_cost, stockout_cost)
  results_ss[[sim]] <- simulate_sq_policy(horizon, s_2, S_2, mean_daily_demand, sd_daily_demand,
                                          lead_time, ordering_cost, holding_cost, stockout_cost)
}

# Aggregate results
sq_results <- tibble(
  total_cost = sapply(results_sq, function(x) x$total_cost),
  service_level = sapply(results_sq, function(x) x$service_level),
  max_inventory = sapply(results_sq, function(x) x$max_inventory)
) |>
  mutate(policy = "(s, Q)")

ss_results <- tibble(
  total_cost = sapply(results_ss, function(x) x$total_cost),
  service_level = sapply(results_ss, function(x) x$service_level),
  max_inventory = sapply(results_ss, function(x) x$max_inventory)
) |>
  mutate(policy = "(s, S)")

combined_results <- bind_rows(sq_results, ss_results)

# Summary statistics
summary_stats <- combined_results |>
  group_by(policy) |>
  summarise(
    mean_cost = mean(total_cost),
    sd_cost = sd(total_cost),
    mean_service_level = 100 * mean(service_level),
    mean_max_inv = mean(max_inventory),
    .groups = "drop"
  )

cat("=== POLICY COMPARISON ===\n\n")
cat("Policy Parameters:\n")
cat(sprintf("(s, Q): s=%d, Q=%d\n", s_1, Q_1))
cat(sprintf("(s, S): s=%d, S=%d\n\n", s_2, S_2))

cat("Performance Metrics (averaged across simulations):\n\n")
print(summary_stats |>
  mutate(
    mean_cost = format(round(mean_cost), big.mark=","),
    sd_cost = format(round(sd_cost), big.mark=","),
    mean_max_inv = format(round(mean_max_inv), big.mark=",")
  ))

cat("\n\nCost Advantage of (s, S) over (s, Q):\n")
sq_cost <- filter(summary_stats, policy == "(s, Q)") |> pull(mean_cost)
ss_cost <- filter(summary_stats, policy == "(s, S)") |> pull(mean_cost)
cost_diff <- sq_cost - ss_cost

if (cost_diff > 0) {
  cat(sprintf("(s, S) saves ₦%.0f per year (%.1f%% lower)\n",
      abs(cost_diff), 100 * abs(cost_diff) / sq_cost))
} else {
  cat(sprintf("(s, Q) saves ₦%.0f per year (%.1f%% lower)\n",
      abs(cost_diff), 100 * abs(cost_diff) / ss_cost))
}

# Distribution plots
library(ggplot2)

ggplot(combined_results, aes(x = total_cost, fill = policy)) +
  geom_histogram(alpha = 0.6, bins = 30) +
  labs(title = "Distribution of Total Inventory Cost",
       x = "Annual Cost (₦)", y = "Frequency") +
  theme_minimal()

ggplot(combined_results, aes(x = service_level, fill = policy)) +
  geom_histogram(alpha = 0.6, bins = 30) +
  labs(title = "Distribution of Service Level",
       x = "Service Level", y = "Frequency") +
  theme_minimal()

ggplot(combined_results, aes(x = policy, y = total_cost, fill = policy)) +
  geom_boxplot(alpha = 0.7) +
  geom_point(alpha = 0.1) +
  labs(title = "Total Cost Comparison",
       x = "Policy", y = "Annual Cost (₦)") +
  theme_minimal() +
  theme(legend.position = "none")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(8293)

# Simulation parameters
n_simulations = 1000
horizon = 365  # days

mean_daily_demand = 50
sd_daily_demand = 8
lead_time = 14
ordering_cost = 100000
holding_cost = 50
stockout_cost = 500

# Policy parameters
s_1, Q_1 = 500, 600   # (s, Q) policy
s_2, S_2 = 500, 1200  # (s, S) policy

def simulate_inventory_policy(n_days, s, Q, mean_d, sd_d, L, order_cost,
                              hold_cost, stockout_cost, is_ss=False, S=None):
    """Simulate inventory policy (s, Q) or (s, S)"""
    inventory = Q
    on_order = 0
    days_until_arrival = 0

    total_holding_cost = 0
    total_ordering_cost = 0
    total_stockout_cost = 0
    num_stockouts = 0
    max_inventory = 0

    for day in range(n_days):
        # Check if order arrives
        if on_order > 0 and days_until_arrival == 0:
            inventory += on_order
            on_order = 0
        elif on_order > 0:
            days_until_arrival -= 1

        # Daily demand
        demand = max(0, np.random.normal(mean_d, sd_d))

        # Fulfill demand
        if inventory >= demand:
            inventory -= demand
        else:
            shortage = demand - inventory
            total_stockout_cost += shortage * stockout_cost
            num_stockouts += 1
            inventory = 0

        # Holding cost
        total_holding_cost += (inventory * hold_cost / 365)
        max_inventory = max(max_inventory, inventory)

        # Reorder
        if inventory < s and on_order == 0:
            if is_ss:
                # (s, S) policy: order up to S
                on_order = S - inventory
            else:
                # (s, Q) policy: order Q
                on_order = Q

            days_until_arrival = L
            total_ordering_cost += order_cost

    service_level = 1 - (num_stockouts / n_days)
    total_cost = total_holding_cost + total_ordering_cost + total_stockout_cost

    return {
        'total_cost': total_cost,
        'service_level': service_level,
        'max_inventory': max_inventory
    }

# Run simulations
print("=== INVENTORY POLICY SIMULATION ===")
print(f"Simulations: {n_simulations} runs × {horizon} days\n")

results_sq = []
results_ss = []

for sim in range(n_simulations):
    results_sq.append(simulate_inventory_policy(
        horizon, s_1, Q_1, mean_daily_demand, sd_daily_demand, lead_time,
        ordering_cost, holding_cost, stockout_cost, is_ss=False
    ))
    results_ss.append(simulate_inventory_policy(
        horizon, s_2, S_2, mean_daily_demand, sd_daily_demand, lead_time,
        ordering_cost, holding_cost, stockout_cost, is_ss=True, S=S_2
    ))

# Convert to DataFrames
sq_df = pd.DataFrame(results_sq)
sq_df['policy'] = '(s, Q)'

ss_df = pd.DataFrame(results_ss)
ss_df['policy'] = '(s, S)'

combined_df = pd.concat([sq_df, ss_df], ignore_index=True)

# Summary statistics
summary = combined_df.groupby('policy').agg({
    'total_cost': ['mean', 'std'],
    'service_level': 'mean',
    'max_inventory': 'mean'
}).round(2)

print("=== POLICY COMPARISON ===\n")
print(f"Policy Parameters:")
print(f"(s, Q): s={s_1}, Q={Q_1}")
print(f"(s, S): s={s_2}, S={S_2}\n")
print("Performance Metrics:\n")
print(summary)

sq_cost = sq_df['total_cost'].mean()
ss_cost = ss_df['total_cost'].mean()
cost_diff = sq_cost - ss_cost

print(f"\n\nCost Comparison:")
if cost_diff > 0:
    print(f"(s, S) saves ₦{abs(cost_diff):,.0f} per year ({100*abs(cost_diff)/sq_cost:.1f}% lower)")
else:
    print(f"(s, Q) saves ₦{abs(cost_diff):,.0f} per year ({100*abs(cost_diff)/ss_cost:.1f}% lower)")

# Visualizations
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Cost distribution
ax = axes[0, 0]
sq_df['total_cost'].hist(bins=30, alpha=0.6, label='(s, Q)', ax=ax, color='blue')
ss_df['total_cost'].hist(bins=30, alpha=0.6, label='(s, S)', ax=ax, color='orange')
ax.set_xlabel('Annual Cost (₦)')
ax.set_ylabel('Frequency')
ax.set_title('Total Cost Distribution')
ax.legend()

# Service level distribution
ax = axes[0, 1]
sq_df['service_level'].hist(bins=30, alpha=0.6, label='(s, Q)', ax=ax, color='blue')
ss_df['service_level'].hist(bins=30, alpha=0.6, label='(s, S)', ax=ax, color='orange')
ax.set_xlabel('Service Level')
ax.set_ylabel('Frequency')
ax.set_title('Service Level Distribution')
ax.legend()

# Boxplot: cost
ax = axes[1, 0]
combined_df.boxplot(column='total_cost', by='policy', ax=ax)
ax.set_xlabel('Policy')
ax.set_ylabel('Annual Cost (₦)')
ax.set_title('Total Cost Comparison')
plt.sca(ax)
plt.xticks(rotation=0)

# Boxplot: service level
ax = axes[1, 1]
combined_df.boxplot(column='service_level', by='policy', ax=ax)
ax.set_xlabel('Policy')
ax.set_ylabel('Service Level')
ax.set_title('Service Level Comparison')

plt.tight_layout()
plt.show()
Caution📝 Section 49.6 Review Questions

  1. What is Monte Carlo simulation? Why is simulation preferable to closed-form analysis for complex inventory policies?

  2. Explain the difference between (s, Q) and (s, S) policies. When would (s, S) lead to lower maximum inventory than (s, Q)?

  3. In the simulation code, what does days_until_arrival track? Why is lead time modeled this way?

  4. A simulation shows that policy A has lower average cost but higher variance in cost. Which policy would a risk-averse manager prefer? Why?

  5. The simulation assumes demand is normally distributed. What if demand is highly skewed (e.g., exponential)? Would the policy parameters change? How would you verify this?

54.7 Case Study: Safety Stock Optimisation for Nigerian Pharmaceutical Distributor

Background: A wholesaler distributes 500 SKUs to 100+ pharmacies and clinics across Nigeria. Inventory value: ₦500m. Carrying cost: ₦125m/year (25%). Stockouts cause customer loss (estimated ₦20/unit margin loss) and goodwill damage (₦500 per stockout event).

Data: 2-year historical sales by SKU; supplier lead times; customer demand patterns by geography (Lagos demand more stable, rural demand erratic).

Analysis:

54.7.1 Case Study Part 1: Synthetic Data Generation and ABC-XYZ

library(tidyverse)

# Generate synthetic case study data for top 50 A-items
set.seed(6541)

pharma_case_data <- tibble(
  sku_id = paste0("A", 1:50),
  sku_name = paste0("Product_", 1:50),
  daily_demand_mean = rnorm(50, mean = 100, sd = 40) |> pmax(10),
  daily_demand_sd = daily_demand_mean * runif(50, min = 0.15, max = 0.4),
  unit_cost = rlnorm(50, meanlog = log(1500), sdlog = 0.6),
  lead_time_days = sample(c(14, 21, 28), size = 50, replace = TRUE),
  supplier = sample(c("supplier_A", "supplier_B", "supplier_C"), size = 50, replace = TRUE)
) |>
  mutate(
    annual_demand = daily_demand_mean * 365,
    annual_value = annual_demand * unit_cost,
    demand_cv = daily_demand_sd / daily_demand_mean,
    xyz_class = case_when(
      demand_cv < 0.5 ~ "X",
      demand_cv < 1.0 ~ "Y",
      TRUE ~ "Z"
    )
  )

cat("=== CASE STUDY: Nigerian Pharmaceutical Distributor ===\n\n")
cat(sprintf("Top 50 A-Items Summary:\n"))
cat(sprintf("Total Annual Demand: %.0f units\n", sum(pharma_case_data$annual_demand)))
cat(sprintf("Total Inventory Value: ₦%.0fm\n\n", sum(pharma_case_data$annual_value) / 1e6))

# Summary by XYZ class
xyz_summary <- pharma_case_data |>
  group_by(xyz_class) |>
  summarise(
    n_items = n(),
    total_value_m = sum(annual_value) / 1e6,
    avg_cv = mean(demand_cv),
    .groups = "drop"
  )

print(xyz_summary)

# Display sample items
cat("\n\nSample A-Items:\n")
print(pharma_case_data |> slice(1:10) |>
  select(sku_id, sku_name, daily_demand_mean, demand_cv, xyz_class, unit_cost, lead_time_days))
import pandas as pd
import numpy as np

np.random.seed(6541)

# Generate synthetic case study data
pharma_case_data = pd.DataFrame({
    'sku_id': [f'A{i}' for i in range(1, 51)],
    'sku_name': [f'Product_{i}' for i in range(1, 51)],
    'daily_demand_mean': np.random.normal(100, 40, 50).clip(10),
    'unit_cost': np.random.lognormal(np.log(1500), 0.6, 50),
    'lead_time_days': np.random.choice([14, 21, 28], size=50),
})

pharma_case_data['daily_demand_sd'] = pharma_case_data['daily_demand_mean'] * np.random.uniform(0.15, 0.4, 50)
pharma_case_data['annual_demand'] = pharma_case_data['daily_demand_mean'] * 365
pharma_case_data['annual_value'] = pharma_case_data['annual_demand'] * pharma_case_data['unit_cost']
pharma_case_data['demand_cv'] = pharma_case_data['daily_demand_sd'] / pharma_case_data['daily_demand_mean']

pharma_case_data['xyz_class'] = pd.cut(pharma_case_data['demand_cv'],
                                        bins=[0, 0.5, 1.0, np.inf],
                                        labels=['X', 'Y', 'Z'],
                                        include_lowest=True)

print("=== CASE STUDY: Nigerian Pharmaceutical Distributor ===\n")
print(f"Top 50 A-Items Summary:")
print(f"Total Annual Demand: {pharma_case_data['annual_demand'].sum():.0f} units")
print(f"Total Inventory Value: ₦{pharma_case_data['annual_value'].sum()/1e6:.0f}m\n")

# Summary by XYZ class
xyz_summary = pharma_case_data.groupby('xyz_class').agg({
    'sku_id': 'count',
    'annual_value': 'sum',
    'demand_cv': 'mean'
}).rename(columns={'sku_id': 'n_items', 'annual_value': 'total_value'})
xyz_summary['total_value_m'] = xyz_summary['total_value'] / 1e6

print("Summary by XYZ Class:")
print(xyz_summary[['n_items', 'total_value_m', 'demand_cv']])

# Display sample items
print("\n\nSample A-Items:")
print(pharma_case_data[['sku_id', 'sku_name', 'daily_demand_mean', 'demand_cv',
                         'xyz_class', 'unit_cost', 'lead_time_days']].head(10).to_string(index=False))

54.7.2 Case Study Part 2: EOQ and Safety Stock for Top Items

# Continue from pharma_case_data

# EOQ and safety stock calculation for all 50 items
ordering_cost <- 100000  # ₦ per order
holding_cost_pct <- 0.25 # 25% of unit cost annually

pharma_analysis <- pharma_case_data |>
  mutate(
    holding_cost = unit_cost * holding_cost_pct,
    eoq = sqrt(2 * annual_demand * ordering_cost / holding_cost),
    orders_per_year = annual_demand / eoq,
    avg_inventory = eoq / 2,

    # Safety stock for 95% service level (z = 1.645)
    z_95 = 1.645,
    sigma_lt_demand_95 = sqrt((daily_demand_sd^2 * lead_time_days) + (daily_demand_mean^2 * (1.5^2))),
    safety_stock_95 = z_95 * sigma_lt_demand_95,
    reorder_point_95 = (daily_demand_mean * lead_time_days) + safety_stock_95,

    # Safety stock for 99% service level (z = 2.326)
    z_99 = 2.326,
    safety_stock_99 = z_99 * sigma_lt_demand_95,
    reorder_point_99 = (daily_demand_mean * lead_time_days) + safety_stock_99
  )

# Top 10 A-items analysis
cat("=== EOQ AND SAFETY STOCK FOR TOP 10 A-ITEMS ===\n\n")

top_10 <- pharma_analysis |>
  arrange(desc(annual_value)) |>
  slice(1:10)

print(top_10 |>
  select(sku_id, annual_demand, unit_cost, eoq, orders_per_year, safety_stock_95, reorder_point_95) |>
  mutate(across(where(is.numeric), round, 0)))

# Scatter plots
library(ggplot2)

ggplot(pharma_analysis, aes(x = unit_cost, y = eoq, size = annual_demand, colour = xyz_class)) +
  geom_point(alpha = 0.6) +
  scale_size_continuous(range = c(2, 8)) +
  scale_colour_manual(values = c("X" = "green", "Y" = "orange", "Z" = "red")) +
  labs(title = "EOQ vs Unit Cost (bubble = annual demand)",
       x = "Unit Cost (₦)", y = "EOQ (units)", colour = "Variability", size = "Annual Demand") +
  theme_minimal()

ggplot(pharma_analysis, aes(x = demand_cv, y = safety_stock_95, colour = xyz_class)) +
  geom_point(size = 3, alpha = 0.7) +
  scale_colour_manual(values = c("X" = "green", "Y" = "orange", "Z" = "red")) +
  labs(title = "Safety Stock vs Demand Variability (95% SL)",
       x = "Coefficient of Variation", y = "Safety Stock (units)", colour = "Class") +
  theme_minimal()

# Total inventory value analysis
total_eoq_value <- sum(pharma_analysis$eoq * pharma_analysis$unit_cost)
total_ss_value_95 <- sum(pharma_analysis$safety_stock_95 * pharma_analysis$unit_cost)
total_ss_value_99 <- sum(pharma_analysis$safety_stock_99 * pharma_analysis$unit_cost)

cat("\n=== TOTAL INVENTORY VALUE (Top 50 A-Items) ===\n")
cat(sprintf("Base Inventory (EOQ/2): ₦%.0fm\n", (sum(pharma_analysis$avg_inventory * pharma_analysis$unit_cost)) / 1e6))
cat(sprintf("Safety Stock (95% SL): ₦%.0fm\n", total_ss_value_95 / 1e6))
cat(sprintf("Safety Stock (99% SL): ₦%.0fm\n", total_ss_value_99 / 1e6))
cat(sprintf("\nAnnual Holding Cost Impact:\n"))
cat(sprintf("  95% SL: ₦%.1fm\n", (total_ss_value_95 * 0.25) / 1e6))
cat(sprintf("  99% SL: ₦%.1fm\n", (total_ss_value_99 * 0.25) / 1e6))
cat(sprintf("  Extra cost for 4% SL increase: ₦%.1fm\n",
    ((total_ss_value_99 - total_ss_value_95) * 0.25) / 1e6))
import numpy as np
from scipy.stats import norm

# EOQ and safety stock calculation
ordering_cost = 100000
holding_cost_pct = 0.25

pharma_analysis = pharma_case_data.copy()
pharma_analysis['holding_cost'] = pharma_analysis['unit_cost'] * holding_cost_pct
pharma_analysis['eoq'] = np.sqrt(2 * pharma_analysis['annual_demand'] * ordering_cost /
                                 pharma_analysis['holding_cost'])
pharma_analysis['orders_per_year'] = pharma_analysis['annual_demand'] / pharma_analysis['eoq']
pharma_analysis['avg_inventory'] = pharma_analysis['eoq'] / 2

# Safety stock for different service levels
lead_time_sd = 1.5  # days
pharma_analysis['sigma_lt_demand'] = np.sqrt(
    (pharma_analysis['daily_demand_sd']**2 * pharma_analysis['lead_time_days']) +
    (pharma_analysis['daily_demand_mean']**2 * lead_time_sd**2)
)

for sl, z in [(95, 1.645), (99, 2.326)]:
    pharma_analysis[f'safety_stock_{sl}'] = z * pharma_analysis['sigma_lt_demand']
    pharma_analysis[f'reorder_point_{sl}'] = (pharma_analysis['daily_demand_mean'] *
                                               pharma_analysis['lead_time_days']) + pharma_analysis[f'safety_stock_{sl}']

# Top 10 A-items
print("=== EOQ AND SAFETY STOCK FOR TOP 10 A-ITEMS ===\n")
top_10 = pharma_analysis.nlargest(10, 'annual_value')
print(top_10[['sku_id', 'annual_demand', 'unit_cost', 'eoq', 'orders_per_year',
               'safety_stock_95', 'reorder_point_95']].round(0).to_string(index=False))

# Visualizations
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# EOQ vs Unit Cost
colors = {'X': 'green', 'Y': 'orange', 'Z': 'red'}
for xyz in ['X', 'Y', 'Z']:
    mask = pharma_analysis['xyz_class'] == xyz
    ax1.scatter(pharma_analysis.loc[mask, 'unit_cost'],
               pharma_analysis.loc[mask, 'eoq'],
               s=pharma_analysis.loc[mask, 'annual_demand']/10,
               alpha=0.6, label=xyz, color=colors[xyz])

ax1.set_xlabel('Unit Cost (₦)')
ax1.set_ylabel('EOQ (units)')
ax1.set_title('EOQ vs Unit Cost')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Safety Stock vs Demand Variability
for xyz in ['X', 'Y', 'Z']:
    mask = pharma_analysis['xyz_class'] == xyz
    ax2.scatter(pharma_analysis.loc[mask, 'demand_cv'],
               pharma_analysis.loc[mask, 'safety_stock_95'],
               s=100, alpha=0.7, label=xyz, color=colors[xyz])

ax2.set_xlabel('Coefficient of Variation')
ax2.set_ylabel('Safety Stock (units)')
ax2.set_title('Safety Stock vs Demand Variability (95% SL)')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Total inventory value analysis
total_eoq_value = (pharma_analysis['eoq']/2 * pharma_analysis['unit_cost']).sum()
total_ss_value_95 = (pharma_analysis['safety_stock_95'] * pharma_analysis['unit_cost']).sum()
total_ss_value_99 = (pharma_analysis['safety_stock_99'] * pharma_analysis['unit_cost']).sum()

print("\n=== TOTAL INVENTORY VALUE (Top 50 A-Items) ===")
print(f"Base Inventory (EOQ/2): ₦{total_eoq_value/1e6:.0f}m")
print(f"Safety Stock (95% SL): ₦{total_ss_value_95/1e6:.0f}m")
print(f"Safety Stock (99% SL): ₦{total_ss_value_99/1e6:.0f}m")
print(f"\nAnnual Holding Cost Impact:")
print(f"  95% SL: ₦{(total_ss_value_95 * 0.25)/1e6:.1f}m")
print(f"  99% SL: ₦{(total_ss_value_99 * 0.25)/1e6:.1f}m")
print(f"  Extra cost for 4% SL increase: ₦{((total_ss_value_99 - total_ss_value_95) * 0.25)/1e6:.1f}m")
Caution📝 Section 49.7 Review Questions

  1. Using the synthetic case study data, identify which A-items are classified as AZ (high value, erratic demand). What inventory management strategy would you recommend for these items?

  2. Calculate the total annual holding cost for safety stock at 95% and 99% service levels for the top 50 A-items. What is the cost of improving from 95% to 99% SL?

  3. A manager proposes reducing the safety stock for all C-items to zero (relying only on mean demand × lead time for reorder point). What are the risks and benefits?

  4. The case study shows that many A-items are supplied by multiple vendors with different lead times. How would you adapt the LP transportation model to account for this?

  5. Propose three data collection improvements that would refine the case study analysis (e.g., actual vs forecast demand, supplier reliability, customer priorities).

54.8 Further Reading

  • Donald Waters (2003). “Inventory Control and Management”
  • Nahmias & Olsen (2015). “Production and Operations Analysis” (7th ed.). Chapters 4–6 cover EOQ, safety stock, and inventory systems.
  • PuLP Library (Python): Linear programming for Python at https://coin-or.github.io/pulp/

Chapter 49 Exercises

1. [Recall] Inventory Costs in Context

Define holding cost, ordering cost, and stockout cost. Give a specific Nigerian business example of each (e.g., a pharmaceutical distributor, beverage wholesaler, or electronics retailer). Estimate the annual cost of each for a distributor holding ₦500m inventory with 20% carrying cost, ₦100k per order, and ₦50/unit stockout loss.

2. [Recall] ABC-XYZ Classification

Explain the ABC classification (by value) and XYZ classification (by demand variability). Create a hypothetical 3×3 ABC-XYZ matrix for a 200-SKU retailer. In which cell would you place a high-volume, low-margin item (e.g., rice)? In which cell a low-volume, high-margin item (e.g., specialty pharmaceuticals)?

3. [Recall] EOQ Formula and Assumptions

State the EOQ formula \(\text{EOQ} = \sqrt{\frac{2Ds}{h}}\) and list its four key assumptions. For each assumption, give a real-world scenario where it is violated in Nigeria (e.g., foreign exchange volatility, port delays, insecurity disrupting supply routes).

4. [Application] EOQ Calculation and Total Cost

A food distributor has annual demand D = 8,000 units, ordering cost s = ₦80,000, and holding cost h = ₦40/unit-year. Calculate: (a) EOQ, (b) number of orders per year, (c) average inventory, (d) total annual cost. If the distributor currently orders 500 units at a time, how much extra does this cost versus EOQ?

5. [Application] Safety Stock and Reorder Point

A pharmacy stocks an antibiotic with daily demand mean = 25 units, SD = 4 units. Lead time is mean = 8 days, SD = 1.5 days. Calculate safety stock and reorder point for (a) 95% service level and (b) 99% service level. If holding cost is ₦20/unit-year, what is the annual extra cost of 99% vs 95% SL for this item alone?

6. [Application] ABC-XYZ Segmentation with Real Data

Using the synthetic dataset generated in Section 49.2 code (or your own 100-SKU dataset), perform ABC-XYZ classification. Create a summary table showing: number of SKUs, total annual value, and average CV for each of the 9 segments. Which segment(s) warrant the highest inventory management effort?

7. [Analysis] Service Level Tradeoff Analysis

A company’s top A-item has annual demand = 10,000 units, unit cost = ₦500, holding cost = 25% annually. Ordering cost = ₦150,000. Stockout cost is estimated at ₦100 per unit of shortage. Calculate EOQ and total cost. Then calculate total annual holding cost for safety stock at service levels of 90%, 95%, and 99%. At what service level is the marginal cost of safety stock approximately equal to the marginal stockout savings?

8. [Analysis] LP Transportation Problem with Capacity

Extend the transportation problem in Section 49.5 to include warehouse capacity constraints: Lagos can hold max ₦250m, Abuja max ₦150m, Kano max ₦100m (valued at selling price). How would you model this in an LP? What new constraints would you add? Solve the modified problem and compare to the unconstrained solution.

9. [Extension] Comparing (s, Q) and (s, S) Policies

Using the simulation code in Section 49.6, conduct a sensitivity analysis: vary the reorder point s from 300 to 700 (in steps of 100) for both (s, Q) and (s, S) policies. Plot average cost, service level, and max inventory as a function of s. Under what conditions does (s, S) outperform (s, Q) on all three metrics simultaneously?

10. [Extension] Demand Forecasting for AZ Items

AZ items (high value, erratic demand) are difficult to manage with standard EOQ-based policies. Propose a machine learning approach to forecast demand for AZ items: (a) describe two suitable models (e.g., time series, regression), (b) discuss validation metrics appropriate for sporadic/intermittent demand, (c) explain how forecast accuracy would feed into safety stock calculation, (d) outline a deployment strategy (quarterly reforecasting, trigger-based updates).

54.9 Chapter 49 Appendix: Mathematical Derivations

54.9.1 A49.1 EOQ Derivation from Calculus

Total cost: \(C(Q) = \frac{Q}{2} h + \frac{D}{Q} s\)

First-order condition: \(\frac{dC}{dQ} = \frac{h}{2} - \frac{Ds}{Q^2} = 0\)

Solving: \(Q^* = \sqrt{\frac{2Ds}{h}}\) (Economic Order Quantity)

Second derivative: \(\frac{d^2C}{dQ^2} = \frac{2Ds}{Q^3} > 0\), confirming minimum.

54.9.2 A49.2 Safety Stock from Normal Distribution

Demand during lead time: \(\mathcal{D}_{LT} \sim \mathcal{N}(\mu_{LT}, \sigma_{LT}^2)\)

\(\mu_{LT} = \mu_d \times L\), \(\sigma_{LT}^2 = \sigma_d^2 \times L + \mu_d^2 \times \sigma_L^2\) (assuming independent demand and lead time)

Service level α: \(P(\text{Inventory} \geq \mathcal{D}_{LT}) = \alpha\), implying \(\text{Inventory} = \mu_{LT} + z_\alpha \sigma_{LT}\)

Safety stock: \(SS = z_\alpha \sigma_{LT}\)

54.9.3 A49.3 ABC-XYZ Segmentation Using Coefficient of Variation

Demand coefficient of variation: \(CV = \frac{\sigma_d}{\mu_d}\)

  • X: CV < 0.5 (stable)
  • Y: 0.5 ≤ CV < 1.0 (variable)
  • Z: CV ≥ 1.0 (erratic)

Pareto A-B-C by value: cumulative revenue % determines thresholds (e.g., top 20% of SKUs = A).

54.9.4 A49.4 Transportation Problem as Linear Program

Variables: \(x_{ij}\) = units from warehouse \(i\) to customer \(j\)

Objective: Minimise \(\sum_i \sum_j c_{ij} x_{ij}\)

Constraints: - Supply: \(\sum_j x_{ij} = S_i\) ∀ \(i\) - Demand: \(\sum_i x_{ij} = D_j\) ∀ \(j\) - Non-negativity: \(x_{ij} ≥ 0\)

The constraint matrix has special structure (network), enabling efficient solutions (network simplex, cost-scaling algorithms).

54.9.5 A49.5 Monte Carlo Simulation for Inventory Policies

For inventory policy (s, Q) over horizon T:

for simulation run = 1 to N:
    inventory ← Q  // start with one order
    total_cost ← 0
    stockouts ← 0
    for t = 1 to T:
        demand_t ← sample from demand distribution
        if inventory < demand_t:
            stockouts ← stockouts + 1
            shortage ← demand_t - inventory
            cost ← shortage × stockout_cost
        else:
            cost ← 0
        inventory ← max(0, inventory - demand_t)
        total_cost ← total_cost + inventory × holding_cost + cost
        if inventory < s:
            inventory ← inventory + Q  // place order
            // advance time by lead_time
    record (total_cost, stockouts, max_inventory)
average across N runs

Output: mean cost, service level (P(no stockout)), max inventory needed.

End of Chapter 49