Supply Chain Simulation in Python — Deep Dive

Production supply chain simulation combines discrete-event modeling with statistical experiment design to quantify risk, size buffers, and validate network redesigns before capital is committed. This guide builds a full multi-tier simulation and walks through disruption modeling, experiment orchestration, and output analysis.

Multi-tier supply chain model

import simpy
import random
import numpy as np
from dataclasses import dataclass, field

@dataclass
class InventoryNode:
    name: str
    env: simpy.Environment
    initial_stock: int
    reorder_point: int
    order_quantity: int
    lead_time_mean: float
    lead_time_std: float
    stock: int = 0
    backlog: int = 0
    total_demand: int = 0
    total_fulfilled: int = 0
    history: list = field(default_factory=list)
    supplier: "InventoryNode | None" = None

    def __post_init__(self):
        self.stock = self.initial_stock

    def fulfill(self, qty: int) -> int:
        served = min(qty, self.stock)
        self.stock -= served
        self.backlog += qty - served
        self.total_demand += qty
        self.total_fulfilled += served
        return served

    def receive(self, qty: int):
        # First clear backlog, then add to stock
        backlog_fill = min(qty, self.backlog)
        self.backlog -= backlog_fill
        self.total_fulfilled += backlog_fill
        self.stock += qty - backlog_fill

    def replenishment_process(self):
        while True:
            yield self.env.timeout(1)  # check daily
            self.history.append({"day": self.env.now, "stock": self.stock, "backlog": self.backlog})

            if self.stock <= self.reorder_point and self.supplier:
                order_qty = self.order_quantity
                lead_time = max(1, random.gauss(self.lead_time_mean, self.lead_time_std))
                self.env.process(self._await_delivery(order_qty, lead_time))

    def _await_delivery(self, qty: int, lead_time: float):
        yield self.env.timeout(lead_time)
        if self.supplier:
            actually_shipped = self.supplier.fulfill(qty)
            self.receive(actually_shipped)
        else:
            self.receive(qty)  # unlimited raw material source

    @property
    def service_level(self) -> float:
        return self.total_fulfilled / max(1, self.total_demand)

Demand generator with seasonality

def demand_generator(env: simpy.Environment, retailer: InventoryNode, base_rate: float, seasonality_amplitude: float = 0.3):
    """Generate daily demand with weekly seasonality and random noise."""
    while True:
        day_of_week = int(env.now) % 7
        seasonal_factor = 1 + seasonality_amplitude * np.sin(2 * np.pi * day_of_week / 7)
        daily_demand = max(0, int(np.random.poisson(base_rate * seasonal_factor)))
        retailer.fulfill(daily_demand)
        yield env.timeout(1)

Disruption modeling

Real supply chains face disruptions: factory fires, port closures, supplier bankruptcies. Model these as stochastic events:

@dataclass
class Disruption:
    name: str
    probability_per_day: float
    duration_mean: float
    duration_std: float
    affected_node: InventoryNode
    capacity_during: float = 0.0  # 0 = complete shutdown

def disruption_process(env: simpy.Environment, disruption: Disruption):
    while True:
        # Time until next disruption (geometric distribution)
        if random.random() < disruption.probability_per_day:
            duration = max(1, random.gauss(disruption.duration_mean, disruption.duration_std))
            original_order_qty = disruption.affected_node.order_quantity
            disruption.affected_node.order_quantity = int(original_order_qty * disruption.capacity_during)
            yield env.timeout(duration)
            disruption.affected_node.order_quantity = original_order_qty
        else:
            yield env.timeout(1)

Example: a port closure with 0.5% daily probability, lasting 7-21 days, that blocks all inbound shipments to a regional warehouse.

Assembling and running the simulation

def run_supply_chain_sim(
    days: int = 365,
    base_demand: float = 50,
    seed: int = 42,
) -> dict:
    random.seed(seed)
    np.random.seed(seed)
    env = simpy.Environment()

    # Three-tier chain: supplier -> warehouse -> retailer
    supplier = InventoryNode(
        name="Supplier", env=env, initial_stock=5000,
        reorder_point=1000, order_quantity=2000,
        lead_time_mean=14, lead_time_std=3, supplier=None,
    )
    warehouse = InventoryNode(
        name="Warehouse", env=env, initial_stock=1000,
        reorder_point=300, order_quantity=500,
        lead_time_mean=5, lead_time_std=1, supplier=supplier,
    )
    retailer = InventoryNode(
        name="Retailer", env=env, initial_stock=200,
        reorder_point=80, order_quantity=150,
        lead_time_mean=2, lead_time_std=0.5, supplier=warehouse,
    )

    # Start processes
    env.process(supplier.replenishment_process())
    env.process(warehouse.replenishment_process())
    env.process(retailer.replenishment_process())
    env.process(demand_generator(env, retailer, base_demand))

    # Add disruption
    port_closure = Disruption(
        name="Port Closure", probability_per_day=0.005,
        duration_mean=14, duration_std=5,
        affected_node=warehouse, capacity_during=0.0,
    )
    env.process(disruption_process(env, port_closure))

    env.run(until=days)

    return {
        "retailer_service_level": retailer.service_level,
        "warehouse_service_level": warehouse.service_level,
        "retailer_avg_stock": np.mean([h["stock"] for h in retailer.history]),
        "retailer_max_backlog": max(h["backlog"] for h in retailer.history) if retailer.history else 0,
    }

Experiment orchestration

Run hundreds of replications with different seeds and parameter variations:

from concurrent.futures import ProcessPoolExecutor
import pandas as pd

def run_experiment(params: dict) -> dict:
    result = run_supply_chain_sim(**params)
    result.update(params)
    return result

experiments = []
for safety_stock_multiplier in [0.5, 1.0, 1.5, 2.0]:
    for seed in range(100):
        experiments.append({
            "days": 365,
            "base_demand": 50,
            "seed": seed,
        })

with ProcessPoolExecutor(max_workers=8) as pool:
    results = list(pool.map(run_experiment, experiments))

df = pd.DataFrame(results)

Analyzing simulation output

Key analyses:

Service level distribution

import matplotlib.pyplot as plt

df.groupby("safety_stock_multiplier")["retailer_service_level"].hist(
    bins=20, alpha=0.5, legend=True
)
plt.xlabel("Service Level")
plt.title("Service Level Distribution Across Scenarios")

Value at Risk (VaR) for stockouts

stockout_days = df.groupby("seed")["retailer_max_backlog"].max()
var_95 = stockout_days.quantile(0.95)
print(f"95% VaR: max backlog of {var_95:.0f} units")

Sensitivity analysis

Which parameter has the most impact on service level? Use correlation or a simple linear regression across experiment results:

from scipy.stats import spearmanr

for param in ["reorder_point", "order_quantity", "lead_time_mean"]:
    corr, p = spearmanr(df[param], df["retailer_service_level"])
    print(f"{param}: Spearman r={corr:.3f}, p={p:.4f}")

Advanced techniques

Agent-based modeling

When supply chain actors make independent decisions (retailers choose suppliers based on price, suppliers allocate scarce stock based on customer priority), agent-based modeling captures emergent behavior. Python’s Mesa library supports agent-based simulation with visualization.

Digital twin integration

Connect the simulation to real-time data feeds (ERP inventory levels, GPS truck positions, weather APIs). The simulation runs continuously, comparing predicted versus actual states and triggering alerts when divergence exceeds thresholds.

Reinforcement learning for policy optimization

Instead of manually setting reorder points and order quantities, train an RL agent to make inventory decisions within the simulation. Libraries like Stable Baselines3 can learn policies that adapt to changing conditions better than static rules.

Pitfalls

• Over-parameterized models — adding every possible disruption type without data to calibrate probabilities produces fiction, not insight. Start with 2-3 well-calibrated disruptions.
• Ignoring warm-up periods — the first 30-60 simulation days reflect initial conditions, not steady state. Discard them from analysis.
• Too few replications — 10 runs are not enough to estimate tail probabilities. Use at least 100 replications; for rare events (1% disruption probability), use 1,000+.
• Confusing the model with reality — a simulation is a simplification. Validate against historical data before trusting its predictions for novel scenarios.

The one thing to remember: Production supply chain simulation chains SimPy process models with disruption events, runs hundreds of Monte Carlo replications, and analyzes the output distribution to quantify risk and guide decisions that spreadsheet models cannot capture.