Stochastic Inventory Models
You are an expert in stochastic inventory theory and probabilistic inventory optimization. Your goal is to help model and optimize inventory systems under demand uncertainty, determining optimal policies that balance inventory costs with service level requirements.
Initial Assessment
Before modeling stochastic inventory, understand:
- Demand Uncertainty
- Demand distribution? (normal, Poisson, negative binomial, empirical) - Demand parameters (mean, variance, coefficient of variation)? - Time period for demand (daily, weekly)? - Intermittent or smooth demand pattern? - Historical data available?
- Lead Time
- Lead time from order to receipt? - Lead time variability? - Lead time distribution? - Correlation between demand and lead time?
- Inventory Policy Type
- Continuous vs. periodic review? - (Q,r) continuous review policy? - (R,S) periodic review policy? - (s,S) policy with bandwidth? - Base stock policy?
- Service Level Requirements
- Target service level? (Type I or Type II) - Type I: Probability of not stocking out during lead time - Type II: Fill rate (fraction of demand satisfied) - Critical vs. non-critical items?
- Cost Structure
- Fixed ordering cost? - Holding cost per unit per period? - Backorder cost vs. lost sales? - Emergency replenishment options?
Stochastic Inventory Fundamentals
Key Concepts
Demand During Lead Time (DDLT):
- Random variable representing total demand during replenishment lead time
- Critical for determining reorder point and safety stock
Safety Stock (SS):
- Buffer inventory to protect against demand uncertainty
- SS = k × σ_DDLT, where k is safety factor
Service Levels:
- Type I (Cycle Service Level, CSL): P(no stockout during lead time)
- Type II (Fill Rate, FR): Fraction of demand met from stock
Inventory Position:
- On-hand + on-order - backorders
- Decision based on position, not just on-hand
Python Implementation: Stochastic Models
(Q,r) Continuous Review Policy
import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import minimize_scalar, fsolve
from typing import Dict, Tuple
import matplotlib.pyplot as plt
class ContinuousReviewQrPolicy:
"""
(Q,r) Continuous Review Inventory Policy
- Monitor inventory continuously
- When inventory position ≤ r, order Q units
- Optimize Q and r to minimize costs while meeting service level
"""
def __init__(self, demand_mean: float, demand_std: float,
lead_time: float, holding_cost: float,
ordering_cost: float, backorder_cost: float = None,
service_level: float = 0.95):
"""
Parameters:
-----------
demand_mean : float
Mean demand per period
demand_std : float
Standard deviation of demand per period
lead_time : float
Lead time in periods
holding_cost : float
Holding cost per unit per period
ordering_cost : float
Fixed ordering cost
backorder_cost : float, optional
Backorder cost per unit per period (if None, uses lost sales)
service_level : float
Target Type I service level (cycle service level)
"""
self.mu = demand_mean
self.sigma = demand_std
self.L = lead_time
self.h = holding_cost
self.K = ordering_cost
self.b = backorder_cost
self.alpha = service_level
# Lead time demand distribution
self.mu_L = demand_mean * lead_time
self.sigma_L = demand_std * np.sqrt(lead_time)
def calculate_eoq(self) -> float:
"""Calculate basic EOQ (ignoring uncertainty)"""
return np.sqrt(2 * self.mu * self.K / self.h)
def calculate_reorder_point(self, Q: float, service_level: float = None) -> Dict:
"""
Calculate reorder point for given order quantity
Parameters:
-----------
Q : float
Order quantity
service_level : float, optional
Target service level (uses self.alpha if None)
Returns:
--------
Dictionary with reorder point, safety stock, and metrics
"""
if service_level is None:
service_level = self.alpha
# Safety factor for target service level
z = stats.norm.ppf(service_level)
# Safety stock
safety_stock = z * self.sigma_L
# Reorder point
r = self.mu_L + safety_stock
# Expected shortage per cycle (for Type II service level)
def G(k):
"""Unit normal loss function"""
return stats.norm.pdf(k) - k * (1 - stats.norm.cdf(k))
k = (r - self.mu_L) / self.sigma_L
expected_shortage = self.sigma_L * G(k)
# Fill rate (Type II service level)
fill_rate = 1 - expected_shortage / Q
return {
'reorder_point': r,
'safety_stock': safety_stock,
'safety_factor': z,
'expected_shortage_per_cycle': expected_shortage,
'fill_rate': fill_rate,
'cycle_service_level': service_level
}
def optimize_Qr_cost(self) -> Dict:
"""
Optimize Q and r jointly to minimize expected total cost
Uses iterative approach:
1. Start with EOQ
2. Calculate optimal r given Q
3. Calculate optimal Q given r
4. Iterate until convergence
"""
# Initialize with EOQ and basic r
Q = self.calculate_eoq()
# Iterate
for _ in range(10):
# Calculate r given Q
r_result = self.calculate_reorder_point(Q, self.alpha)
r = r_result['reorder_point']
# Calculate Q given r (considering safety stock)
# Modified EOQ accounting for safety stock
Q_new = np.sqrt(2 * self.mu * self.K / self.h)
if abs(Q_new - Q) < 0.01:
break
Q = Q_new
# Final calculations
r_result = self.calculate_reorder_point(Q, self.alpha)
r = r_result['reorder_point']
ss = r_result['safety_stock']
# Cost calculations
ordering_cost_annual = (self.mu / Q) * self.K
holding_cost_annual = (Q / 2 + ss) * self.h
# Backorder cost (if specified)
if self.b is not None:
expected_shortage = r_result['expected_shortage_per_cycle']
backorder_cost_annual = (self.mu / Q) * expected_shortage * self.b
total_cost = ordering_cost_annual + holding_cost_annual + backorder_cost_annual
else:
backorder_cost_annual = 0
total_cost = ordering_cost_annual + holding_cost_annual
return {
'policy': '(Q,r)',
'order_quantity': Q,
'reorder_point': r,
'safety_stock': ss,
'safety_factor': r_result['safety_factor'],
'avg_inventory': Q / 2 + ss,
'ordering_cost': ordering_cost_annual,
'holding_cost': holding_cost_annual,
'backorder_cost': backorder_cost_annual,
'total_cost': total_cost,
'cycle_service_level': r_result['cycle_service_level'],
'fill_rate': r_result['fill_rate']
}
def simulate(self, num_periods: int, Q: float, r: float,
seed: int = None) -> pd.DataFrame:
"""
Simulate (Q,r) policy over time
Parameters:
-----------
num_periods : int
Number of periods to simulate
Q : float
Order quantity
r : float
Reorder point
seed : int, optional
Random seed for reproducibility
"""
if seed is not None:
np.random.seed(seed)
# Initialize
inventory_position = Q # Start with one order
inventory_on_hand = Q
orders_outstanding = [] # List of (arrival_period, quantity)
results = []
for t in range(num_periods):
# Receive orders arriving this period
arrivals = [order for order in orders_outstanding
if order['arrival'] == t]
for arrival in arrivals:
inventory_on_hand += arrival['quantity']
orders_outstanding.remove(arrival)
# Generate demand
demand = max(0, np.random.normal(self.mu, self.sigma))
# Satisfy demand
sales = min(demand, inventory_on_hand)
stockout = demand - sales
inventory_on_hand -= sales
# Update inventory position
inventory_position = inventory_on_hand + sum(
o['quantity'] for o in orders_outstanding)
# Check if order needed
order_placed = 0
if inventory_position <= r:
order_placed = Q
orders_outstanding.append({
'order_period': t,
'arrival': t + int(self.L),
'quantity': Q
})
inventory_position += Q
results.append({
'period': t,
'demand': demand,
'sales': sales,
'stockout': stockout,
'inventory_on_hand': inventory_on_hand,
'inventory_position': inventory_position,
'order_placed': order_placed,
'orders_outstanding': len(orders_outstanding)
})
df = pd.DataFrame(results)
# Calculate performance metrics
fill_rate = df['sales'].sum() / df['demand'].sum()
avg_inventory = df['inventory_on_hand'].mean()
stockout_periods = (df['stockout'] > 0).sum()
cycle_service_level = 1 - stockout_periods / num_periods
summary = {
'fill_rate': fill_rate,
'cycle_service_level': cycle_service_level,
'avg_inventory': avg_inventory,
'total_stockouts': df['stockout'].sum(),
'num_orders': (df['order_placed'] > 0).sum()
}
return df, summary
def plot_simulation(self, simulation_df: pd.DataFrame, Q: float, r: float):
"""Visualize simulation results"""
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
periods = simulation_df['period']
# Plot 1: Inventory levels and reorder point
axes[0].plot(periods, simulation_df['inventory_position'],
label='Inventory Position', linewidth=2, color='blue')
axes[0].plot(periods, simulation_df['inventory_on_hand'],
label='On-Hand Inventory', linewidth=2, color='green', alpha=0.7)
axes[0].axhline(y=r, color='red', linestyle='--', linewidth=2,
label=f'Reorder Point (r={r:.0f})')
axes[0].fill_between(periods, 0, simulation_df['inventory_on_hand'],
alpha=0.2, color='green')
# Mark order placements
order_periods = periods[simulation_df['order_placed'] > 0]
axes[0].scatter(order_periods,
simulation_df.loc[simulation_df['order_placed'] > 0,
'inventory_position'],
color='red', s=100, marker='^', zorder=5,
label='Order Placed')
axes[0].set_ylabel('Inventory Level')
axes[0].set_title('(Q,r) Policy: Inventory Levels Over Time', fontweight='bold')
axes[0].legend(loc='upper right')
axes[0].grid(True, alpha=0.3)
# Plot 2: Demand and sales
axes[1].plot(periods, simulation_df['demand'], label='Demand',
linewidth=2, color='blue', alpha=0.6)
axes[1].plot(periods, simulation_df['sales'], label='Sales',
linewidth=2, color='green')
axes[1].fill_between(periods, simulation_df['sales'],
simulation_df['demand'],
where=(simulation_df['stockout'] > 0),
alpha=0.3, color='red', label='Stockouts')
axes[1].set_ylabel('Units')
axes[1].set_title('Demand vs. Sales', fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
# Plot 3: Cumulative stockouts
axes[2].plot(periods, simulation_df['stockout'].cumsum(),
linewidth=2, color='red')
axes[2].set_xlabel('Period')
axes[2].set_ylabel('Cumulative Stockout Units')
axes[2].set_title('Cumulative Stockouts Over Time', fontweight='bold')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
return plt
# Example Usage
def example_continuous_review():
"""Example: (Q,r) continuous review policy"""
print("\n" + "=" * 70)
print("(Q,r) CONTINUOUS REVIEW STOCHASTIC INVENTORY MODEL")
print("=" * 70)
# Create model
model = ContinuousReviewQrPolicy(
demand_mean=100, # Average demand: 100 units/day
demand_std=20, # Std dev: 20 units/day
lead_time=7, # Lead time: 7 days
holding_cost=2.0, # $2/unit/day
ordering_cost=200, # $200 per order
backorder_cost=10, # $10/unit backordered
service_level=0.95 # 95% cycle service level
)
print("\nInput Parameters:")
print(f" Demand: Normal(μ={model.mu}, σ={model.sigma}) per day")
print(f" Lead Time: {model.L} days")
print(f" Lead Time Demand: Normal(μ={model.mu_L:.0f}, σ={model.sigma_L:.1f})")
print(f" Holding Cost: ${model.h}/unit/day")
print(f" Ordering Cost: ${model.K}")
print(f" Backorder Cost: ${model.b}/unit/day")
print(f" Target Service Level: {model.alpha:.1%}")
# Optimize policy
optimal = model.optimize_Qr_cost()
print(f"\n{'=' * 70}")
print("OPTIMAL (Q,r) POLICY")
print("=" * 70)
print(f"\n{'Order Quantity (Q):':<35} {optimal['order_quantity']:.0f} units")
print(f"{'Reorder Point (r):':<35} {optimal['reorder_point']:.0f} units")
print(f"{'Safety Stock:':<35} {optimal['safety_stock']:.0f} units")
print(f"{'Safety Factor (z):':<35} {optimal['safety_factor']:.2f}")
print(f"{'Average Inventory:':<35} {optimal['avg_inventory']:.0f} units")
print(f"\n{'Performance Metrics:':<35}")
print(f" {'Cycle Service Level:':<33} {optimal['cycle_service_level']:.1%}")
print(f" {'Fill Rate:':<33} {optimal['fill_rate']:.2%}")
print(f"\n{'Annual Costs:':<35}")
print(f" {'Ordering Cost:':<33} ${optimal['ordering_cost']:,.2f}")
print(f" {'Holding Cost:':<33} ${optimal['holding_cost']:,.2f}")
print(f" {'Backorder Cost:':<33} ${optimal['backorder_cost']:,.2f}")
print(f" {'-'*50}")
print(f" {'Total Annual Cost:':<33} ${optimal['total_cost']:,.2f}")
# Simulate
print("\n\nSimulating policy for 365 days...")
simulation, summary = model.simulate(
num_periods=365,
Q=optimal['order_quantity'],
r=optimal['reorder_point'],
seed=42
)
print(f"\n{'Simulation Results (365 days):':<35}")
print(f" {'Fill Rate (actual):':<33} {summary['fill_rate']:.2%}")
print(f" {'Cycle Service Level (actual):':<33} {summary['cycle_service_level']:.1%}")
print(f" {'Average Inventory:':<33} {summary['avg_inventory']:.0f} units")
print(f" {'Total Stockout Units:':<33} {summary['total_stockouts']:.0f}")
print(f" {'Number of Orders:':<33} {summary['num_orders']}")
# Plot
model.plot_simulation(simulation, optimal['order_quantity'],
optimal['reorder_point'])
plt.savefig('/tmp/qr_policy_simulation.png', dpi=300, bbox_inches='tight')
print(f"\nSimulation plot saved to /tmp/qr_policy_simulation.png")
return model, optimal, simulation
if __name__ == "__main__":
example_continuous_review()Periodic Review (R,S) Policy
class PeriodicReviewRSPolicy:
"""
(R,S) Periodic Review Policy
- Review inventory every R periods
- Order up to level S (order-up-to level)
- Must cover demand during R + L periods
"""
def __init__(self, demand_mean: float, demand_std: float,
lead_time: float, review_period: float,
holding_cost: float, ordering_cost: float,
service_level: float = 0.95):
"""
Parameters:
-----------
demand_mean : float
Mean demand per period
demand_std : float
Std dev of demand per period
lead_time : float
Lead time in periods
review_period : float
Time between reviews (R)
holding_cost : float
Holding cost per unit per period
ordering_cost : float
Fixed ordering cost per order
service_level : float
Target Type I service level
"""
self.mu = demand_mean
self.sigma = demand_std
self.L = lead_time
self.R = review_period
self.h = holding_cost
self.K = ordering_cost
self.alpha = service_level
# Protection period = R + L
self.T = review_period + lead_time
# Demand distribution during protection period
self.mu_T = demand_mean * self.T
self.sigma_T = demand_std * np.sqrt(self.T)
def calculate_order_up_to_level(self) -> Dict:
"""
Calculate order-up-to level S
S = E[Demand during R+L] + Safety Stock
"""
# Safety factor
z = stats.norm.ppf(self.alpha)
# Safety stock
safety_stock = z * self.sigma_T
# Order-up-to level
S = self.mu_T + safety_stock
# Average inventory (approximation)
# On average, order R*mu units every R periods
avg_order = self.R * self.mu
avg_inventory = avg_order / 2 + safety_stock
# Costs
orders_per_year = 1 / self.R # If period = year/365
ordering_cost_annual = orders_per_year * self.K
holding_cost_annual = avg_inventory * self.h * self.R
total_cost = ordering_cost_annual + holding_cost_annual
return {
'policy': '(R,S)',
'review_period': self.R,
'order_up_to_level': S,
'safety_stock': safety_stock,
'safety_factor': z,
'avg_inventory': avg_inventory,
'ordering_cost': ordering_cost_annual,
'holding_cost': holding_cost_annual,
'total_cost': total_cost,
'service_level': self.alpha
}
# Example: Periodic Review
def example_periodic_review():
"""Example: (R,S) periodic review policy"""
print("\n" + "=" * 70)
print("(R,S) PERIODIC REVIEW STOCHASTIC INVENTORY MODEL")
print("=" * 70)
model = PeriodicReviewRSPolicy(
demand_mean=50, # 50 units/day
demand_std=10, # Std dev 10
lead_time=5, # 5 days lead time
review_period=7, # Review weekly
holding_cost=1.5, # $1.50/unit/day
ordering_cost=150, # $150 per order
service_level=0.98 # 98% service level
)
print("\nInput Parameters:")
print(f" Review Period (R): {model.R} days")
print(f" Lead Time (L): {model.L} days")
print(f" Protection Period (R+L): {model.T} days")
print(f" Demand per Day: Normal(μ={model.mu}, σ={model.sigma})")
print(f" Demand during Protection: Normal(μ={model.mu_T:.0f}, σ={model.sigma_T:.1f})")
result = model.calculate_order_up_to_level()
print(f"\n{'=' * 70}")
print("OPTIMAL (R,S) POLICY")
print("=" * 70)
print(f"\n{'Review Period (R):':<35} {result['review_period']:.0f} days")
print(f"{'Order-Up-To Level (S):':<35} {result['order_up_to_level']:.0f} units")
print(f"{'Safety Stock:':<35} {result['safety_stock']:.0f} units")
print(f"{'Safety Factor:':<35} {result['safety_factor']:.2f}")
print(f"{'Target Service Level:':<35} {result['service_level']:.1%}")
print(f"\n{'Policy:':<35} Review inventory every {result['review_period']:.0f} days")
print(f"{' ':<35} Order up to {result['order_up_to_level']:.0f} units")
return model, result
if __name__ == "__main__":
example_periodic_review()Tools & Libraries
Python Libraries
scipy.stats: Probability distributionsnumpy,pandas: Numerical and data operations- Custom implementations for policy optimization
Commercial Software
- Blue Yonder: Advanced stochastic inventory optimization
- ToolsGroup: Probabilistic demand forecasting and inventory
- Logility: Inventory optimization with uncertainty
- SAP IBP: Stochastic planning capabilities
- o9 Solutions: Probabilistic inventory planning
Common Challenges & Solutions
Challenge: Intermittent Demand
Problem: Many zeros, high variability Solutions:
- Use specialized distributions (Croston's method, Poisson, negative binomial)
- Higher safety stocks
- Consider make-to-order
Challenge: Lead Time Variability
Problem: Supplier lead times uncertain Solutions:
- Model lead time as random variable
- Combined uncertainty formula: σ²_DDLT = L·σ²_D + μ²_D·σ²_L
- Safety lead time approach
Challenge: Service Level Selection
Problem: Difficult to specify target Solutions:
- Use cost-based optimization (balance holding vs. stockout)
- ABC classification with differentiated service
- Empirically validate with business
Challenge: Demand Distribution Selection
Problem: Don't know appropriate distribution Solutions:
- Use empirical bootstrap methods
- Fit and test multiple distributions (normal, lognormal, gamma)
- Normal often works well for aggregate demand
Related Skills
- inventory-optimization: General inventory management
- economic-order-quantity: Deterministic models
- newsvendor-problem: Single-period stochastic
- multi-echelon-inventory: Network-wide stochastic models
- demand-forecasting: Demand distribution estimation
- safety-stock: Safety stock calculation methods