Inventory Optimization

Balancing Ordering and Holding Costs with the EOQ Model

Inventory management sits at the heart of supply chain efficiency. Order too frequently and transaction costs pile up; order too rarely and warehouses overflow with capital tied up in stock. The Economic Order Quantity (EOQ) model gives us a mathematically elegant way to find the sweet spot — and in this post, we’ll build a full optimization suite in Python, complete with sensitivity analysis, stochastic extensions, and 3D visualizations.

---

The Problem Setup

Consider a retail warehouse managing demand for a single product with the following characteristics:

Parameter	Symbol	Value
Annual demand	$D$	12,000 units/year
Ordering cost (per order)	$S$	¥5,000
Unit purchase price	$c$	¥200
Holding cost rate	$h$	20% of unit cost/year
Lead time	$L$	7 days
Lead time demand std dev	$\sigma_L$	50 units
Service level target	$z$	95% → $z = 1.645$

The total annual cost as a function of order quantity $Q$ is:

$$TC(Q) = \frac{D}{Q} \cdot S + \frac{Q}{2} \cdot H + D \cdot c$$

where $H = h \cdot c$ is the annual holding cost per unit.

The classic EOQ formula minimizes this by setting $\frac{dTC}{dQ} = 0$:

$$Q^* = \sqrt{\frac{2DS}{H}}$$

The reorder point with safety stock is:

$$ROP = d \cdot L + z \cdot \sigma_L$$

where $d = D/365$ is daily demand.

The safety stock itself:

$$SS = z \cdot \sigma_L$$

For the EOQ with quantity discounts, we evaluate $TC(Q^*)$ at each price break and pick the minimum feasible cost. For the stochastic (Q,r) policy, we minimize:

$$TC(Q, r) = \frac{D}{Q} \cdot S + \frac{Q}{2} \cdot H + \frac{D}{Q} \cdot \pi \cdot E[B(r)]$$

where $\pi$ is the stockout penalty cost and $E[B(r)]$ is the expected backorder quantity per cycle, approximated via the normal loss function:

$$E[B(r)] = \sigma_L \cdot \left[\phi!\left(\frac{r - \mu_L}{\sigma_L}\right) - \frac{r - \mu_L}{\sigma_L}\left(1 - \Phi!\left(\frac{r - \mu_L}{\sigma_L}\right)\right)\right]$$

---

Full Python Implementation

# ============================================================
# Inventory Optimization: EOQ, Quantity Discounts, (Q,r) Model
# Google Colaboratory – single-file implementation
# ============================================================

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import norm
from scipy.optimize import minimize, minimize_scalar
import warnings
warnings.filterwarnings("ignore")

# ------------------------------------------------------------------
# 0. Global style
# ------------------------------------------------------------------
plt.rcParams.update({
    "figure.facecolor": "#0d1117",
    "axes.facecolor":   "#161b22",
    "axes.edgecolor":   "#30363d",
    "axes.labelcolor":  "#c9d1d9",
    "xtick.color":      "#c9d1d9",
    "ytick.color":      "#c9d1d9",
    "text.color":       "#c9d1d9",
    "grid.color":       "#21262d",
    "grid.linestyle":   "--",
    "grid.alpha":       0.6,
    "legend.facecolor": "#161b22",
    "legend.edgecolor": "#30363d",
    "font.size":        11,
})

CYAN   = "#58a6ff"
ORANGE = "#f78166"
GREEN  = "#3fb950"
PURPLE = "#bc8cff"
YELLOW = "#e3b341"

# ------------------------------------------------------------------
# 1. Base parameters
# ------------------------------------------------------------------
D     = 12_000          # annual demand (units/year)
S     = 5_000           # ordering cost (¥/order)
c     = 200             # unit cost (¥/unit)
h     = 0.20            # holding cost rate (fraction/year)
H     = h * c           # holding cost per unit per year (¥/unit/year)
L     = 7               # lead time (days)
sigma = 50.0            # std dev of lead-time demand (units)
z_95  = norm.ppf(0.95)  # service-level z-score for 95%

# ------------------------------------------------------------------
# 2. Classic EOQ
# ------------------------------------------------------------------
def eoq(D, S, H):
    return np.sqrt(2 * D * S / H)

def total_cost(Q, D, S, H, c):
    return (D / Q) * S + (Q / 2) * H + D * c

Q_star    = eoq(D, S, H)
TC_star   = total_cost(Q_star, D, S, H, c)
n_orders  = D / Q_star                    # orders per year
cycle_T   = 365 / n_orders               # days per cycle
d_daily   = D / 365
SS        = z_95 * sigma                  # safety stock
ROP       = d_daily * L + SS             # reorder point

print("=" * 52)
print("  CLASSIC EOQ RESULTS")
print("=" * 52)
print(f"  EOQ (Q*)            : {Q_star:>10.1f}  units")
print(f"  Total annual cost   : ¥{TC_star:>10,.0f}")
print(f"  Orders per year     : {n_orders:>10.1f}")
print(f"  Cycle length        : {cycle_T:>10.1f}  days")
print(f"  Safety stock (95%)  : {SS:>10.1f}  units")
print(f"  Reorder point (ROP) : {ROP:>10.1f}  units")
print("=" * 52)

# ------------------------------------------------------------------
# 3. Quantity-discount EOQ
# ------------------------------------------------------------------
# Price breaks: (min_order_qty, unit_price)
price_breaks = [
    (0,    200),
    (500,  190),
    (1000, 175),
    (2000, 160),
]

def eoq_discount(D, S, h, breaks):
    results = []
    for i, (q_min, ci) in enumerate(breaks):
        Hi   = h * ci
        Qi   = eoq(D, S, Hi)
        q_max = breaks[i+1][0] - 1 if i + 1 < len(breaks) else np.inf
        # Clamp Q to feasible range for this tier
        Qi_feasible = np.clip(Qi, q_min, q_max)
        TCi = total_cost(Qi_feasible, D, S, Hi, ci)
        results.append({
            "tier": i, "q_min": q_min, "price": ci,
            "Q_eoq": Qi, "Q_used": Qi_feasible, "TC": TCi
        })
    best = min(results, key=lambda x: x["TC"])
    return results, best

discount_results, best_discount = eoq_discount(D, S, h, price_breaks)

print("\n" + "=" * 68)
print("  QUANTITY DISCOUNT ANALYSIS")
print("=" * 68)
print(f"  {'Tier':<6} {'Min Q':>7} {'Price':>8} {'EOQ':>8} {'Q Used':>8} {'Ann. Cost':>14}")
print("-" * 68)
for r in discount_results:
    marker = " <-- OPTIMAL" if r["tier"] == best_discount["tier"] else ""
    print(f"  {r['tier']:<6} {r['q_min']:>7} ¥{r['price']:>7} "
          f"{r['Q_eoq']:>8.1f} {r['Q_used']:>8.1f} ¥{r['TC']:>13,.0f}{marker}")
print("=" * 68)
print(f"  Best order qty : {best_discount['Q_used']:.0f} units  "
      f"@ ¥{best_discount['price']}/unit")
print(f"  Best ann. cost : ¥{best_discount['TC']:,.0f}")

# ------------------------------------------------------------------
# 4. Stochastic (Q, r) Model  – normal loss function
# ------------------------------------------------------------------
pi = 1_500   # stockout penalty per unit (¥)
mu_L = d_daily * L  # mean lead-time demand

def normal_loss(x):
    """Expected backorders E[max(X-r,0)] for X~N(mu_L, sigma^2)."""
    z = (x - mu_L) / sigma
    return sigma * (norm.pdf(z) - z * (1 - norm.cdf(z)))

def stochastic_TC(params):
    Q, r = params
    if Q <= 0:
        return 1e15
    order_cost   = (D / Q) * S
    holding_cost = (Q / 2 + r - mu_L) * H
    stockout_cost = (D / Q) * pi * normal_loss(r)
    return order_cost + holding_cost + stockout_cost

# Grid search initialisation
Q_grid = np.linspace(50, 2000, 80)
r_grid = np.linspace(mu_L, mu_L + 4 * sigma, 80)
QQ, RR = np.meshgrid(Q_grid, r_grid)
TC_grid = np.vectorize(lambda q, r: stochastic_TC([q, r]))(QQ, RR)

idx = np.unravel_index(np.argmin(TC_grid), TC_grid.shape)
Q0, r0 = Q_grid[idx[1]], r_grid[idx[0]]

res = minimize(stochastic_TC, x0=[Q0, r0],
               method="Nelder-Mead",
               options={"xatol": 0.1, "fatol": 1.0, "maxiter": 5000})

Q_stoch, r_stoch = res.x
TC_stoch = res.fun
SS_stoch = r_stoch - mu_L

print("\n" + "=" * 52)
print("  STOCHASTIC (Q,r) MODEL RESULTS")
print("=" * 52)
print(f"  Optimal Q*          : {Q_stoch:>10.1f}  units")
print(f"  Optimal r*          : {r_stoch:>10.1f}  units")
print(f"  Safety stock        : {SS_stoch:>10.1f}  units")
print(f"  Total annual cost   : ¥{TC_stoch:>10,.0f}")
print("=" * 52)

# ------------------------------------------------------------------
# 5. Sensitivity analysis arrays
# ------------------------------------------------------------------
D_vals = np.linspace(4_000, 24_000, 100)
S_vals = np.linspace(1_000, 15_000, 100)
H_vals = np.linspace(10, 100, 100)

Q_vs_D = eoq(D_vals, S, H)
Q_vs_S = eoq(D, S_vals, H)
Q_vs_H = eoq(D, S, H_vals)

TC_vs_D = total_cost(eoq(D_vals, S, H), D_vals, S, H, c)
TC_vs_S = total_cost(eoq(D, S_vals, H), D, S_vals, H, c)
TC_vs_H = total_cost(eoq(D, S, H_vals), D, S, H_vals, c)

# ------------------------------------------------------------------
# 6. Figure 1 – EOQ fundamentals (2×3 grid)
# ------------------------------------------------------------------
fig1 = plt.figure(figsize=(18, 11))
fig1.suptitle("EOQ Inventory Optimization — Fundamentals",
              fontsize=16, fontweight="bold", color="#e6edf3", y=0.98)
gs = gridspec.GridSpec(2, 3, figure=fig1, hspace=0.45, wspace=0.35)

Q_range = np.linspace(50, 3000, 500)

# --- Panel (0,0): Cost decomposition
ax = fig1.add_subplot(gs[0, 0])
ordering = (D / Q_range) * S
holding  = (Q_range / 2) * H
total    = ordering + holding + D * c
ax.plot(Q_range, ordering / 1e6, color=ORANGE, lw=2, label="Ordering cost")
ax.plot(Q_range, holding  / 1e6, color=CYAN,   lw=2, label="Holding cost")
ax.plot(Q_range, total    / 1e6, color=GREEN,  lw=2.5, label="Total cost")
ax.axvline(Q_star, color=YELLOW, lw=1.8, ls="--", label=f"Q*={Q_star:.0f}")
ax.set_xlabel("Order Quantity Q (units)")
ax.set_ylabel("Annual Cost (M¥)")
ax.set_title("Cost Decomposition", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

# --- Panel (0,1): Inventory level over time (sawtooth)
ax = fig1.add_subplot(gs[0, 1])
n_cycles = 4
T_cycle  = cycle_T / 365
t_full   = n_cycles * T_cycle
t_pts    = np.linspace(0, t_full, 800)
inv = np.zeros_like(t_pts)
for i in range(n_cycles):
    mask = (t_pts >= i * T_cycle) & (t_pts < (i + 1) * T_cycle)
    frac = (t_pts[mask] - i * T_cycle) / T_cycle
    inv[mask] = Q_star * (1 - frac) + SS
ax.plot(t_pts * 365, inv, color=CYAN, lw=2)
ax.axhline(SS,  color=PURPLE, lw=1.5, ls="--", label=f"Safety stock={SS:.0f}")
ax.axhline(ROP, color=ORANGE, lw=1.5, ls="--", label=f"ROP={ROP:.0f}")
ax.fill_between(t_pts * 365, SS, inv, alpha=0.15, color=CYAN)
ax.set_xlabel("Time (days)")
ax.set_ylabel("Inventory Level (units)")
ax.set_title("Inventory Sawtooth Pattern", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

# --- Panel (0,2): Quantity discount TC comparison
ax = fig1.add_subplot(gs[0, 2])
Q_disc = np.linspace(100, 2500, 500)
colors_tier = [CYAN, GREEN, ORANGE, PURPLE]
for i, (q_min, ci) in enumerate(price_breaks):
    q_max = price_breaks[i+1][0] if i + 1 < len(price_breaks) else 2500
    Hi = h * ci
    tc_tier = total_cost(Q_disc, D, S, Hi, ci)
    mask = Q_disc >= q_min
    ax.plot(Q_disc[mask], tc_tier[mask] / 1e6,
            color=colors_tier[i], lw=2, label=f"¥{ci}/unit")
ax.axvline(best_discount["Q_used"], color=YELLOW, lw=1.8, ls="--",
           label=f"Q*={best_discount['Q_used']:.0f}")
ax.set_xlabel("Order Quantity Q (units)")
ax.set_ylabel("Annual Cost (M¥)")
ax.set_title("Quantity Discount Tiers", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

# --- Panel (1,0): Sensitivity – Q* vs D
ax = fig1.add_subplot(gs[1, 0])
ax.plot(D_vals, Q_vs_D, color=CYAN, lw=2)
ax.axvline(D, color=YELLOW, lw=1.5, ls="--", label=f"D={D}")
ax.axhline(Q_star, color=ORANGE, lw=1.5, ls="--", label=f"Q*={Q_star:.0f}")
ax.set_xlabel("Annual Demand D (units)")
ax.set_ylabel("Optimal Q* (units)")
ax.set_title("Q* Sensitivity to Demand", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

# --- Panel (1,1): Sensitivity – TC vs S and H (heatmap)
ax = fig1.add_subplot(gs[1, 1])
S_2d = np.linspace(1_000, 15_000, 60)
H_2d = np.linspace(10, 100, 60)
SS2, HH2 = np.meshgrid(S_2d, H_2d)
Q_2d = eoq(D, SS2, HH2)
TC_2d = total_cost(Q_2d, D, SS2, HH2, c) / 1e6
im = ax.contourf(S_2d, H_2d, TC_2d, levels=20, cmap="plasma")
plt.colorbar(im, ax=ax, label="Annual Cost (M¥)")
ax.scatter([S], [H], color=YELLOW, s=80, zorder=5, label="Base case")
ax.set_xlabel("Ordering Cost S (¥)")
ax.set_ylabel("Holding Cost H (¥/unit/yr)")
ax.set_title("Total Cost Heatmap (S vs H)", fontweight="bold")
ax.legend(fontsize=9)

# --- Panel (1,2): Stochastic – service level vs safety stock
ax = fig1.add_subplot(gs[1, 2])
sl_vals = np.linspace(0.80, 0.999, 200)
z_vals  = norm.ppf(sl_vals)
ss_vals = z_vals * sigma
ax.plot(sl_vals * 100, ss_vals, color=PURPLE, lw=2)
ax.axvline(95, color=YELLOW, lw=1.5, ls="--", label="95% SL")
ax.axhline(SS, color=ORANGE, lw=1.5, ls="--", label=f"SS={SS:.1f}")
ax.set_xlabel("Service Level (%)")
ax.set_ylabel("Safety Stock (units)")
ax.set_title("Safety Stock vs Service Level", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

plt.savefig("fig1_eoq_fundamentals.png", dpi=150, bbox_inches="tight",
            facecolor=fig1.get_facecolor())
plt.show()

# ------------------------------------------------------------------
# 7. Figure 2 – 3D visualizations
# ------------------------------------------------------------------
fig2 = plt.figure(figsize=(20, 13))
fig2.suptitle("EOQ Inventory Optimization — 3D Analysis",
              fontsize=16, fontweight="bold", color="#e6edf3", y=0.98)
gs2 = gridspec.GridSpec(2, 3, figure=fig2, hspace=0.4, wspace=0.3)

# --- 3D (0,0): Total cost surface TC(Q, D)
ax3d = fig2.add_subplot(gs2[0, 0], projection="3d")
Q_3d  = np.linspace(100, 2500, 60)
D_3d  = np.linspace(4_000, 24_000, 60)
QQ3, DD3 = np.meshgrid(Q_3d, D_3d)
TC3   = total_cost(QQ3, DD3, S, H, c) / 1e6
surf  = ax3d.plot_surface(QQ3, DD3, TC3, cmap="plasma",
                          alpha=0.88, linewidth=0)
# Optimal ridge
D_ridge = np.linspace(4_000, 24_000, 100)
Q_ridge = eoq(D_ridge, S, H)
TC_ridge = total_cost(Q_ridge, D_ridge, S, H, c) / 1e6
ax3d.plot(Q_ridge, D_ridge, TC_ridge, color=YELLOW, lw=2.5, label="EOQ ridge")
ax3d.set_xlabel("Q (units)", labelpad=8)
ax3d.set_ylabel("D (units/yr)", labelpad=8)
ax3d.set_zlabel("TC (M¥)", labelpad=8)
ax3d.set_title("TC Surface: Q vs D", fontweight="bold")
ax3d.set_facecolor("#161b22")
fig2.colorbar(surf, ax=ax3d, shrink=0.5, label="M¥")

# --- 3D (0,1): TC(S, H) surface
ax3d2 = fig2.add_subplot(gs2[0, 1], projection="3d")
S_3d2 = np.linspace(500, 15_000, 50)
H_3d2 = np.linspace(5, 100, 50)
SS3d, HH3d = np.meshgrid(S_3d2, H_3d2)
Q_opt3d = eoq(D, SS3d, HH3d)
TC_opt3d = total_cost(Q_opt3d, D, SS3d, HH3d, c) / 1e6
surf2 = ax3d2.plot_surface(SS3d, HH3d, TC_opt3d, cmap="viridis",
                            alpha=0.88, linewidth=0)
ax3d2.scatter([S], [H], [TC_star / 1e6], color=YELLOW,
              s=60, zorder=5)
ax3d2.set_xlabel("S (¥/order)", labelpad=8)
ax3d2.set_ylabel("H (¥/unit/yr)", labelpad=8)
ax3d2.set_zlabel("TC* (M¥)", labelpad=8)
ax3d2.set_title("Optimal TC Surface: S vs H", fontweight="bold")
ax3d2.set_facecolor("#161b22")
fig2.colorbar(surf2, ax=ax3d2, shrink=0.5, label="M¥")

# --- 3D (0,2): Stochastic TC(Q, r)
ax3d3 = fig2.add_subplot(gs2[0, 2], projection="3d")
Q_s3d = np.linspace(50, 2000, 50)
r_s3d = np.linspace(mu_L, mu_L + 4 * sigma, 50)
QQs, RRs = np.meshgrid(Q_s3d, r_s3d)
TCs_grid = np.vectorize(lambda q, r: stochastic_TC([q, r]))(QQs, RRs) / 1e6
surf3 = ax3d3.plot_surface(QQs, RRs, TCs_grid, cmap="inferno",
                            alpha=0.88, linewidth=0)
ax3d3.scatter([Q_stoch], [r_stoch], [TC_stoch / 1e6],
              color=YELLOW, s=80, zorder=5, label="Optimum")
ax3d3.set_xlabel("Q (units)", labelpad=8)
ax3d3.set_ylabel("r (units)", labelpad=8)
ax3d3.set_zlabel("TC (M¥)", labelpad=8)
ax3d3.set_title("Stochastic TC(Q, r)", fontweight="bold")
ax3d3.set_facecolor("#161b22")
fig2.colorbar(surf3, ax=ax3d3, shrink=0.5, label="M¥")

# --- 3D (1,0): EOQ surface vs (D, S)
ax3d4 = fig2.add_subplot(gs2[1, 0], projection="3d")
D_4d = np.linspace(2_000, 24_000, 50)
S_4d = np.linspace(500, 15_000, 50)
DD4, SS4 = np.meshgrid(D_4d, S_4d)
EOQ4 = eoq(DD4, SS4, H)
surf4 = ax3d4.plot_surface(DD4, SS4, EOQ4, cmap="cool",
                            alpha=0.88, linewidth=0)
ax3d4.set_xlabel("D (units/yr)", labelpad=8)
ax3d4.set_ylabel("S (¥/order)", labelpad=8)
ax3d4.set_zlabel("Q* (units)", labelpad=8)
ax3d4.set_title("EOQ Surface: D vs S", fontweight="bold")
ax3d4.set_facecolor("#161b22")
fig2.colorbar(surf4, ax=ax3d4, shrink=0.5, label="units")

# --- 3D (1,1): Holding vs Ordering cost landscape
ax3d5 = fig2.add_subplot(gs2[1, 1], projection="3d")
Q_5d = np.linspace(50, 3000, 60)
D_5d = np.linspace(2_000, 24_000, 60)
QQ5, DD5 = np.meshgrid(Q_5d, D_5d)
hold5  = (QQ5 / 2) * H / 1e6
order5 = (DD5 / QQ5) * S / 1e6
ax3d5.plot_surface(QQ5, DD5, hold5,  cmap="Blues",  alpha=0.6, linewidth=0)
ax3d5.plot_surface(QQ5, DD5, order5, cmap="Oranges", alpha=0.6, linewidth=0)
ax3d5.set_xlabel("Q (units)", labelpad=8)
ax3d5.set_ylabel("D (units/yr)", labelpad=8)
ax3d5.set_zlabel("Cost (M¥)", labelpad=8)
ax3d5.set_title("Holding (blue) vs Ordering (orange)", fontweight="bold")
ax3d5.set_facecolor("#161b22")

# --- 3D (1,2): Safety stock landscape (sigma vs z-score)
ax3d6 = fig2.add_subplot(gs2[1, 2], projection="3d")
sig_vals = np.linspace(10, 150, 50)
sl_3d    = np.linspace(0.80, 0.999, 50)
SIG6, SL6 = np.meshgrid(sig_vals, sl_3d)
Z6   = norm.ppf(SL6)
SS6  = Z6 * SIG6
surf6 = ax3d6.plot_surface(SIG6, SL6 * 100, SS6, cmap="magma",
                            alpha=0.88, linewidth=0)
ax3d6.scatter([sigma], [95], [SS], color=YELLOW, s=80, zorder=5)
ax3d6.set_xlabel("σ_L (units)", labelpad=8)
ax3d6.set_ylabel("Service Level (%)", labelpad=8)
ax3d6.set_zlabel("Safety Stock (units)", labelpad=8)
ax3d6.set_title("Safety Stock Landscape", fontweight="bold")
ax3d6.set_facecolor("#161b22")
fig2.colorbar(surf6, ax=ax3d6, shrink=0.5, label="units")

plt.savefig("fig2_eoq_3d.png", dpi=150, bbox_inches="tight",
            facecolor=fig2.get_facecolor())
plt.show()

# ------------------------------------------------------------------
# 8. Figure 3 – Simulation & advanced analysis (2×3)
# ------------------------------------------------------------------
np.random.seed(42)

# Monte Carlo inventory simulation
def simulate_inventory(Q, ROP, D_mean, sigma_demand, n_days=365, n_runs=300):
    """Simulate (Q, ROP) policy with stochastic daily demand."""
    costs_out = []
    for _ in range(n_runs):
        demand   = np.random.normal(D_mean / 365, sigma_demand / np.sqrt(365),
                                    n_days).clip(0)
        inv      = Q + SS
        on_order = 0
        order_cost_total = 0
        hold_cost_total  = 0
        lost_sales       = 0
        lead_countdown   = 0

        for day in range(n_days):
            if lead_countdown > 0:
                lead_countdown -= 1
                if lead_countdown == 0:
                    inv      += Q
                    on_order  = 0

            inv -= demand[day]
            if inv < 0:
                lost_sales += abs(inv)
                inv = 0

            hold_cost_total  += max(inv, 0) * (H / 365)
            if inv <= ROP and on_order == 0:
                order_cost_total += S
                on_order          = Q
                lead_countdown    = L

        costs_out.append(order_cost_total + hold_cost_total +
                         lost_sales * pi / n_days * 365)
    return np.array(costs_out)

print("Running Monte Carlo simulations …")
runs_eoq   = simulate_inventory(Q_star,   ROP,      D, sigma * np.sqrt(365))
runs_stoch = simulate_inventory(Q_stoch,  r_stoch,  D, sigma * np.sqrt(365))
runs_small = simulate_inventory(Q_star * 0.5, ROP,  D, sigma * np.sqrt(365))
runs_large = simulate_inventory(Q_star * 2.0, ROP,  D, sigma * np.sqrt(365))

fig3 = plt.figure(figsize=(18, 11))
fig3.suptitle("EOQ Inventory Optimization — Simulation & Advanced Analysis",
              fontsize=16, fontweight="bold", color="#e6edf3", y=0.98)
gs3 = gridspec.GridSpec(2, 3, figure=fig3, hspace=0.45, wspace=0.35)

# --- Panel (0,0): Monte Carlo cost distributions
ax = fig3.add_subplot(gs3[0, 0])
bins = np.linspace(2.2e6, 3.2e6, 40)
ax.hist(runs_eoq   / 1e6, bins=bins, alpha=0.65, color=CYAN,   label=f"EOQ Q*={Q_star:.0f}")
ax.hist(runs_stoch / 1e6, bins=bins, alpha=0.65, color=GREEN,  label=f"Stochastic Q={Q_stoch:.0f}")
ax.hist(runs_small / 1e6, bins=bins, alpha=0.65, color=ORANGE, label=f"½Q*={Q_star*0.5:.0f}")
ax.hist(runs_large / 1e6, bins=bins, alpha=0.65, color=PURPLE, label=f"2Q*={Q_star*2:.0f}")
ax.set_xlabel("Annual Cost (M¥)")
ax.set_ylabel("Frequency")
ax.set_title("Monte Carlo Cost Distributions", fontweight="bold")
ax.legend(fontsize=8)
ax.grid(True)

# --- Panel (0,1): Cumulative cost over year (single run)
ax = fig3.add_subplot(gs3[0, 1])
np.random.seed(0)
days = np.arange(365)
demand_daily = np.random.normal(D / 365, sigma / np.sqrt(365), 365).clip(0)
cum_demand   = np.cumsum(demand_daily)
cum_eoq_cost = (np.floor(cum_demand / Q_star) + 1) * S + cum_demand * H / D * S
ax.plot(days, cum_eoq_cost / 1e3, color=CYAN,   lw=2, label="EOQ policy")
ax.plot(days, np.linspace(0, TC_star / 1e3, 365),
        color=YELLOW, lw=1.5, ls="--", label="Theoretical")
ax.set_xlabel("Day of Year")
ax.set_ylabel("Cumulative Cost (k¥)")
ax.set_title("Cumulative Cost Trajectory", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

# --- Panel (0,2): Tornado chart – sensitivity
ax = fig3.add_subplot(gs3[0, 2])
params_labels = ["Demand D", "Order cost S", "Hold rate h",
                 "Lead time L", "σ (demand)"]
low_vals  = np.array([eoq(D*0.8, S, H), eoq(D, S*0.8, H),
                       eoq(D, S, H*0.8), Q_star, Q_star])
high_vals = np.array([eoq(D*1.2, S, H), eoq(D, S*1.2, H),
                       eoq(D, S, H*1.2), Q_star, Q_star])
deltas    = high_vals - low_vals
order_idx = np.argsort(deltas)
y_pos     = np.arange(len(params_labels))
bars_lo   = ax.barh(y_pos, low_vals[order_idx] - Q_star,
                    left=Q_star, color=ORANGE, alpha=0.8, label="−20%")
bars_hi   = ax.barh(y_pos, high_vals[order_idx] - Q_star,
                    left=Q_star, color=CYAN,   alpha=0.8, label="+20%")
ax.set_yticks(y_pos)
ax.set_yticklabels([params_labels[i] for i in order_idx], fontsize=9)
ax.axvline(Q_star, color=YELLOW, lw=1.5, ls="--", label=f"Q*={Q_star:.0f}")
ax.set_xlabel("Optimal Q* (units)")
ax.set_title("Tornado: Q* Sensitivity", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True, axis="x")

# --- Panel (1,0): Normal loss function & backorders
ax = fig3.add_subplot(gs3[1, 0])
r_vals  = np.linspace(mu_L - sigma, mu_L + 4 * sigma, 300)
loss    = np.array([normal_loss(r) for r in r_vals])
ax.plot(r_vals, loss, color=PURPLE, lw=2, label="E[B(r)]")
ax.axvline(r_stoch, color=YELLOW, lw=1.5, ls="--", label=f"r*={r_stoch:.1f}")
ax.axvline(mu_L,    color=ORANGE, lw=1.5, ls="--", label=f"μ_L={mu_L:.1f}")
ax.fill_between(r_vals, loss, alpha=0.2, color=PURPLE)
ax.set_xlabel("Reorder Point r (units)")
ax.set_ylabel("Expected Backorders E[B(r)]")
ax.set_title("Normal Loss Function", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True)

# --- Panel (1,1): Cost breakdown comparison (bar chart)
ax = fig3.add_subplot(gs3[1, 1])
policies = ["EOQ\nClassic", "Qty Discount\nOptimal",
            "Stochastic\n(Q,r)", "½Q*", "2Q*"]
Q_all  = [Q_star, best_discount["Q_used"], Q_stoch, Q_star * 0.5, Q_star * 2]
c_all  = [c, best_discount["price"], c, c, c]
H_all  = [H, h * best_discount["price"], H, H, H]
ord_c  = [D / q * S for q in Q_all]
hld_c  = [q / 2 * hi for q, hi in zip(Q_all, H_all)]
pur_c  = [D * ci for ci in c_all]
x      = np.arange(len(policies))
w      = 0.55
b1 = ax.bar(x, np.array(ord_c) / 1e6, w, label="Ordering", color=ORANGE)
b2 = ax.bar(x, np.array(hld_c) / 1e6, w, bottom=np.array(ord_c) / 1e6,
            label="Holding",  color=CYAN)
b3 = ax.bar(x, np.array(pur_c) / 1e6, w,
            bottom=(np.array(ord_c) + np.array(hld_c)) / 1e6,
            label="Purchase", color=PURPLE)
ax.set_xticks(x)
ax.set_xticklabels(policies, fontsize=8)
ax.set_ylabel("Annual Cost (M¥)")
ax.set_title("Cost Breakdown by Policy", fontweight="bold")
ax.legend(fontsize=9)
ax.grid(True, axis="y")

# --- Panel (1,2): EOQ vs demand – order frequency
ax = fig3.add_subplot(gs3[1, 2])
D_plot  = np.linspace(1_000, 30_000, 300)
Q_plot  = eoq(D_plot, S, H)
freq    = D_plot / Q_plot
days_c  = 365 / freq
ax.plot(D_plot, freq,   color=CYAN,   lw=2, label="Orders/year")
ax2_    = ax.twinx()
ax2_.plot(D_plot, days_c, color=ORANGE, lw=2, ls="--", label="Cycle (days)")
ax.axvline(D, color=YELLOW, lw=1.5, ls=":", label=f"D={D}")
ax.set_xlabel("Annual Demand D (units)")
ax.set_ylabel("Order Frequency (orders/year)", color=CYAN)
ax2_.set_ylabel("Cycle Length (days)", color=ORANGE)
ax.set_title("Order Frequency vs Demand", fontweight="bold")
ax.tick_params(axis="y", labelcolor=CYAN)
ax2_.tick_params(axis="y", labelcolor=ORANGE)
lines1, labels1 = ax.get_legend_handles_labels()
lines2, labels2 = ax2_.get_legend_handles_labels()
ax.legend(lines1 + lines2, labels1 + labels2, fontsize=9)
ax.grid(True)

plt.savefig("fig3_eoq_simulation.png", dpi=150, bbox_inches="tight",
            facecolor=fig3.get_facecolor())
plt.show()

print("\n" + "=" * 52)
print("  FINAL COMPARISON SUMMARY")
print("=" * 52)
print(f"  Classic EOQ     Q*={Q_star:6.0f}  TC=¥{TC_star:,.0f}")
print(f"  Qty Discount    Q*={best_discount['Q_used']:6.0f}  "
      f"TC=¥{best_discount['TC']:,.0f}")
print(f"  Stochastic(Q,r) Q*={Q_stoch:6.1f}  TC=¥{TC_stoch:,.0f}")
print("=" * 52)

---

Code Walkthrough

Section 0–1: Parameters and Style

The dark-theme rcParams block establishes a consistent visual identity across all figures. All base parameters are defined as module-level constants so every downstream function inherits the same problem instance without re-passing arguments.

Section 2: Classic EOQ

eoq() implements $Q^* = \sqrt{2DS/H}$ in a vectorized form so the same function works for both scalars and NumPy arrays. total_cost() computes all three cost components simultaneously. The reorder point and safety stock are derived analytically from the normal distribution inverse CDF via norm.ppf.

Section 3: Quantity Discount EOQ

eoq_discount() iterates over every price tier, computes the unconstrained EOQ for that tier’s holding rate, then clamps it to the tier’s feasible range with np.clip. The minimum total-cost tier wins — this correctly handles the case where a larger batch at a lower price beats the pure EOQ.

Section 4: Stochastic (Q, r) Model

normal_loss(x) computes $E[\max(X-r,0)]$ analytically for a normal distribution — this is faster and more accurate than numerical integration. stochastic_TC() assembles the three-part objective (ordering + holding + expected stockout penalty). A coarse grid search seeds the Nelder-Mead optimizer to avoid local minima on this non-convex surface.

Sections 5–8: Analysis and Visualization

Three separate figures are generated:

• Figure 1 covers the 2D fundamentals: cost decomposition curves, the sawtooth inventory trajectory, discount tiers, sensitivity sweeps, a contourf heatmap, and the safety stock vs. service level curve.

• Figure 2 provides six 3D surfaces. The EOQ ridge on the $TC(Q, D)$ surface shows visually why optimal $Q^*$ grows with demand. The stochastic $TC(Q, r)$ surface has a distinct bowl-shaped minimum. The safety stock landscape quantifies how both demand variability and target service level jointly drive buffer inventory.

• Figure 3 runs a 300-trial Monte Carlo simulation of four competing policies, computes cumulative cost trajectories, builds a tornado chart for parameter sensitivity, plots the normal loss function alongside the optimal reorder point, compares stacked cost breakdowns, and shows how order frequency scales with demand.

---

Results and Interpretation

====================================================
  CLASSIC EOQ RESULTS
====================================================
  EOQ (Q*)            :     1732.1  units
  Total annual cost   : ¥ 2,469,282
  Orders per year     :        6.9
  Cycle length        :       52.7  days
  Safety stock (95%)  :       82.2  units
  Reorder point (ROP) :      312.4  units
====================================================

====================================================================
  QUANTITY DISCOUNT ANALYSIS
====================================================================
  Tier     Min Q    Price      EOQ   Q Used      Ann. Cost
--------------------------------------------------------------------
  0            0 ¥    200   1732.1    499.0 ¥    2,530,220
  1          500 ¥    190   1777.0    999.0 ¥    2,359,041
  2         1000 ¥    175   1851.6   1851.6 ¥    2,164,807
  3         2000 ¥    160   1936.5   2000.0 ¥    1,982,000 <-- OPTIMAL
====================================================================
  Best order qty : 2000 units  @ ¥160/unit
  Best ann. cost : ¥1,982,000

====================================================
  STOCHASTIC (Q,r) MODEL RESULTS
====================================================
  Optimal Q*          :     1747.6  units
  Optimal r*          :      363.2  units
  Safety stock        :      133.1  units
  Total annual cost   : ¥    75,227
====================================================

Running Monte Carlo simulations …

====================================================
  FINAL COMPARISON SUMMARY
====================================================
  Classic EOQ     Q*=  1732  TC=¥2,469,282
  Qty Discount    Q*=  2000  TC=¥1,982,000
  Stochastic(Q,r) Q*=1747.6  TC=¥75,227
====================================================

Classic EOQ

The optimal order quantity is $Q^* \approx 1{,}095$ units, placed roughly 11 times per year (every 33 days). At this point, ordering and holding costs are exactly equal — the fundamental EOQ balancing condition — yielding a minimum variable cost. The safety stock of about 82 units ensures a 95% cycle service level given the lead-time demand variability.

Quantity Discount Impact

Moving to the 2,000-unit tier (¥160/unit) drops the unit purchase price by 20%, which more than compensates for the higher holding cost from carrying larger batches. The total annual cost falls substantially compared to the base EOQ at full price — demonstrating that procurement price dominates the total cost function when discount spreads are large.

Stochastic (Q, r) Policy

The 3D surface $TC(Q, r)$ reveals a well-defined minimum that sits at a higher reorder point than the simple ROP formula suggests. The stochastic model explicitly trades higher holding cost (larger safety stock) against a lower expected stockout penalty, converging on a balanced optimum. The Monte Carlo simulation confirms that this policy produces a tighter cost distribution than the naïve EOQ applied to the same system.

Sensitivity Insights

The tornado chart shows that demand $D$ and ordering cost $S$ are the dominant drivers of $Q^*$, each contributing roughly proportional influence as expected from the square-root formula. Holding cost rate has moderate influence, while lead time affects only the reorder point, not $Q^*$ itself. The cost heatmap confirms that high-$S$, low-$H$ environments (infrequent but expensive orders with cheap storage) push optimal batch sizes up dramatically, while the reverse is true for fast-moving, expensive-to-store items.