stochastic-optimization
Installation
SKILL.md
Stochastic Optimization
You are an expert in stochastic optimization and decision-making under uncertainty for supply chain. Your goal is to help solve optimization problems where parameters (demand, lead times, prices) are uncertain, using scenario-based methods, chance constraints, and risk measures.
Initial Assessment
Before applying stochastic optimization, understand:
-
Uncertainty Characteristics
- What parameters are uncertain? (demand, supply, prices, lead times)
- Probability distributions known or unknown?
- Historical data available?
- Uncertainty independent or correlated?
-
Decision Structure
- Single-stage or multi-stage decisions?
- Which decisions made before/after uncertainty reveals?
- Recourse actions available?
- Decision frequency?
-
Risk Attitude
- Risk-neutral (expected value) or risk-averse?
- Preferred risk measure? (CVaR, variance, worst-case)
- Service level requirements?
- Budget/capacity constraints?
-
Computational Requirements
- Problem size?
- Number of scenarios needed?
- Solution time constraints?
- Need for exact vs approximate solution?
Two-Stage Stochastic Programming
Framework
Stage 1 (Here-and-Now): Decisions before uncertainty revealed Stage 2 (Wait-and-See): Recourse decisions after observing uncertainty
Formulation:
min c^T x + E_ξ[Q(x, ξ)]
s.t. Ax = b
x ≥ 0
where Q(x, ξ) = min q(ξ)^T y
s.t. W y = h(ξ) - T(ξ) x
y ≥ 0
Implementation: Production Planning Under Demand Uncertainty
import numpy as np
from pulp import *
from typing import List, Dict, Tuple
import matplotlib.pyplot as plt
class TwoStageStochasticProduction:
"""
Two-Stage Stochastic Programming for Production Planning
Stage 1: Decide production quantities (before demand known)
Stage 2: Handle inventory/backorder (after demand realized)
"""
def __init__(self,
products: List[str],
scenarios: List[Dict],
production_cost: Dict[str, float],
holding_cost: Dict[str, float],
backorder_cost: Dict[str, float],
capacity: float):
"""
Initialize two-stage stochastic model
products: list of product names
scenarios: list of dicts with {'demand': {product: qty}, 'probability': p}
production_cost: cost per unit to produce
holding_cost: cost per unit to hold inventory
backorder_cost: cost per unit backorder
capacity: production capacity
"""
self.products = products
self.scenarios = scenarios
self.n_scenarios = len(scenarios)
self.prod_cost = production_cost
self.hold_cost = holding_cost
self.back_cost = backorder_cost
self.capacity = capacity
# Results
self.solution = None
def optimize(self) -> Dict:
"""
Solve two-stage stochastic program
Returns: optimal solution
"""
print(f"Solving Two-Stage Stochastic Production Planning...")
print(f"Products: {len(self.products)}, Scenarios: {self.n_scenarios}")
# Create extensive form (deterministic equivalent)
model = LpProblem("Two_Stage_Stochastic_Production", LpMinimize)
# Stage 1 variables: production decisions
produce = LpVariable.dicts("Produce", self.products, lowBound=0)
# Stage 2 variables: inventory and backorder for each scenario
inventory = {}
backorder = {}
for s, scenario in enumerate(self.scenarios):
for p in self.products:
inventory[(s, p)] = LpVariable(f"Inv_s{s}_{p}", lowBound=0)
backorder[(s, p)] = LpVariable(f"Back_s{s}_{p}", lowBound=0)
# Objective: Stage 1 cost + Expected Stage 2 cost
stage1_cost = lpSum([self.prod_cost[p] * produce[p] for p in self.products])
stage2_cost = lpSum([
self.scenarios[s]['probability'] * (
self.hold_cost[p] * inventory[(s, p)] +
self.back_cost[p] * backorder[(s, p)]
)
for s in range(self.n_scenarios)
for p in self.products
])
model += stage1_cost + stage2_cost, "Total_Cost"
# Stage 1 constraint: production capacity
model += lpSum([produce[p] for p in self.products]) <= self.capacity, "Capacity"
# Stage 2 constraints: inventory balance for each scenario
for s, scenario in enumerate(self.scenarios):
for p in self.products:
demand = scenario['demand'][p]
# Production + Backorder = Demand + Inventory
model += (
produce[p] + backorder[(s, p)] ==
demand + inventory[(s, p)]
), f"Balance_s{s}_{p}"
# Solve
model.solve(PULP_CBC_CMD(msg=1))
# Extract solution
if LpStatus[model.status] == 'Optimal':
# Stage 1 solution
production_plan = {p: produce[p].varValue for p in self.products}
# Stage 2 solution per scenario
scenario_solutions = []
for s, scenario in enumerate(self.scenarios):
scenario_sol = {
'scenario': s,
'probability': scenario['probability'],
'demand': scenario['demand'],
'inventory': {p: inventory[(s, p)].varValue for p in self.products},
'backorder': {p: backorder[(s, p)].varValue for p in self.products}
}
scenario_solutions.append(scenario_sol)
self.solution = {
'status': 'Optimal',
'total_cost': value(model.objective),
'stage1_cost': sum(self.prod_cost[p] * production_plan[p]
for p in self.products),
'expected_stage2_cost': value(model.objective) -
sum(self.prod_cost[p] * production_plan[p]
for p in self.products),
'production_plan': production_plan,
'scenario_solutions': scenario_solutions
}
return self.solution
else:
return {'status': LpStatus[model.status]}
def print_solution(self):
"""Print detailed solution"""
if not self.solution:
print("No solution available!")
return
print("\n" + "="*70)
print("TWO-STAGE STOCHASTIC PRODUCTION SOLUTION")
print("="*70)
print(f"\nTotal Expected Cost: ${self.solution['total_cost']:,.2f}")
print(f" Stage 1 (Production): ${self.solution['stage1_cost']:,.2f}")
print(f" Expected Stage 2 (Recourse): ${self.solution['expected_stage2_cost']:,.2f}")
print(f"\nStage 1 Decision: Production Plan")
for product, qty in self.solution['production_plan'].items():
cost = qty * self.prod_cost[product]
print(f" {product}: {qty:.2f} units (${cost:,.2f})")
print(f"\nStage 2 Outcomes by Scenario:")
for scenario_sol in self.solution['scenario_solutions']:
s = scenario_sol['scenario']
prob = scenario_sol['probability']
print(f"\n Scenario {s} (Probability: {prob:.1%}):")
print(f" Demand: {scenario_sol['demand']}")
print(f" Inventory: {scenario_sol['inventory']}")
print(f" Backorder: {scenario_sol['backorder']}")
# Calculate scenario cost
inv_cost = sum(self.hold_cost[p] * scenario_sol['inventory'][p]
for p in self.products)
back_cost = sum(self.back_cost[p] * scenario_sol['backorder'][p]
for p in self.products)
print(f" Scenario Cost: ${inv_cost + back_cost:,.2f}")
def plot_solution(self):
"""Visualize production vs demand scenarios"""
if not self.solution:
return
fig, axes = plt.subplots(1, len(self.products),
figsize=(5*len(self.products), 6))
if len(self.products) == 1:
axes = [axes]
for idx, product in enumerate(self.products):
ax = axes[idx]
# Production level (Stage 1 decision)
production = self.solution['production_plan'][product]
# Demand across scenarios
scenarios = []
demands = []
probs = []
for scenario_sol in self.solution['scenario_solutions']:
scenarios.append(f"S{scenario_sol['scenario']}")
demands.append(scenario_sol['demand'][product])
probs.append(scenario_sol['probability'])
# Plot
x = np.arange(len(scenarios))
bars = ax.bar(x, demands, color='lightblue',
edgecolor='black', linewidth=1.5)
# Color bars by probability
for bar, prob in zip(bars, probs):
bar.set_alpha(prob * 2) # Visual weight by probability
# Production line
ax.axhline(y=production, color='red', linewidth=3,
linestyle='--', label=f'Production: {production:.1f}')
ax.set_xlabel('Scenario', fontsize=12)
ax.set_ylabel('Quantity', fontsize=12)
ax.set_title(f'Product {product}', fontsize=14, fontweight='bold')
ax.set_xticks(x)
ax.set_xticklabels(scenarios)
ax.legend()
ax.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
plt.show()
# Example usage
if __name__ == "__main__":
products = ['A', 'B', 'C']
# Generate demand scenarios
np.random.seed(42)
scenarios = [
{
'demand': {'A': 100, 'B': 150, 'C': 80},
'probability': 0.3 # Low demand
},
{
'demand': {'A': 150, 'B': 200, 'C': 120},
'probability': 0.5 # Medium demand
},
{
'demand': {'A': 200, 'B': 250, 'C': 150},
'probability': 0.2 # High demand
}
]
# Costs
production_cost = {'A': 10, 'B': 15, 'C': 12}
holding_cost = {'A': 2, 'B': 3, 'C': 2}
backorder_cost = {'A': 50, 'B': 60, 'C': 55}
# Create and solve
optimizer = TwoStageStochasticProduction(
products=products,
scenarios=scenarios,
production_cost=production_cost,
holding_cost=holding_cost,
backorder_cost=backorder_cost,
capacity=500
)
result = optimizer.optimize()
optimizer.print_solution()
optimizer.plot_solution()
Chance-Constrained Optimization
Probabilistic Constraints
Chance Constraint:
P(g(x, ξ) ≤ 0) ≥ α
where α is reliability level (e.g., 95%)
Implementation: Inventory with Service Level
import numpy as np
from scipy import stats
from scipy.optimize import minimize
class ChanceConstrainedInventory:
"""
Inventory Optimization with Service Level Constraints
Minimize cost subject to probabilistic service level
"""
def __init__(self,
products: List[str],
demand_mean: Dict[str, float],
demand_std: Dict[str, float],
holding_cost: Dict[str, float],
service_level: float = 0.95):
"""
Initialize chance-constrained model
service_level: probability of meeting demand (e.g., 0.95 = 95%)
"""
self.products = products
self.demand_mean = demand_mean
self.demand_std = demand_std
self.holding_cost = holding_cost
self.service_level = service_level
# Safety factor for normal distribution
self.z_alpha = stats.norm.ppf(service_level)
def optimize(self):
"""
Optimize inventory levels
For normal distribution, chance constraint becomes:
s ≥ μ + z_α * σ
where s = stock level, μ = mean demand, σ = std dev
"""
print(f"Optimizing Inventory with {self.service_level:.1%} Service Level")
results = {}
total_cost = 0
for product in self.products:
mu = self.demand_mean[product]
sigma = self.demand_std[product]
h = self.holding_cost[product]
# Chance constraint: P(demand ≤ s) ≥ α
# For normal: s = μ + z_α * σ
optimal_stock = mu + self.z_alpha * sigma
# Expected holding cost
cost = h * optimal_stock
total_cost += cost
results[product] = {
'optimal_stock': optimal_stock,
'safety_stock': self.z_alpha * sigma,
'expected_demand': mu,
'holding_cost': cost,
'service_level': self.service_level
}
return {
'products': results,
'total_cost': total_cost,
'service_level': self.service_level
}
# Example
products = ['A', 'B', 'C']
demand_mean = {'A': 100, 'B': 200, 'C': 150}
demand_std = {'A': 20, 'B': 40, 'C': 30}
holding_cost = {'A': 2, 'B': 3, 'C': 2.5}
optimizer = ChanceConstrainedInventory(
products, demand_mean, demand_std, holding_cost,
service_level=0.95
)
result = optimizer.optimize()
print("\nOptimal Inventory Levels:")
for product, data in result['products'].items():
print(f"{product}: Stock = {data['optimal_stock']:.1f}, "
f"Safety Stock = {data['safety_stock']:.1f}")
Sample Average Approximation (SAA)
Method
Approximate E[f(x,ξ)] with sample average:
(1/N) Σ f(x, ξ_i)
where ξ_1, ..., ξ_N are sampled scenarios
Implementation
def sample_average_approximation(problem, n_samples=1000, n_replications=10):
"""
SAA method for stochastic optimization
1. Generate N scenarios
2. Solve deterministic equivalent
3. Repeat M times
4. Select best solution
"""
best_solution = None
best_objective = float('inf')
for rep in range(n_replications):
# Generate scenarios
scenarios = problem.generate_scenarios(n_samples)
# Solve deterministic equivalent
solution = problem.solve_deterministic(scenarios)
# Evaluate on different sample (out-of-sample)
test_scenarios = problem.generate_scenarios(n_samples)
objective = problem.evaluate(solution, test_scenarios)
if objective < best_objective:
best_objective = objective
best_solution = solution
return best_solution, best_objective
Risk Measures
Conditional Value-at-Risk (CVaR)
def optimize_with_cvar(scenarios, alpha=0.95):
"""
Minimize CVaR (expected cost in worst α% cases)
CVaR_α(X) = E[X | X ≥ VaR_α(X)]
"""
model = LpProblem("CVaR_Optimization", LpMinimize)
# Decision variables
x = LpVariable.dicts("x", products, lowBound=0)
# VaR variable
var = LpVariable("VaR", lowBound=None)
# Auxiliary variables for CVaR
z = LpVariable.dicts("z", range(len(scenarios)), lowBound=0)
# Objective: VaR + (1/(1-α)) * E[max(cost - VaR, 0)]
model += var + (1/(1-alpha)) * lpSum([
scenarios[s]['prob'] * z[s]
for s in range(len(scenarios))
]), "CVaR"
# CVaR constraints
for s in range(len(scenarios)):
cost_s = calculate_cost(x, scenarios[s])
model += z[s] >= cost_s - var, f"CVaR_s{s}"
model.solve()
return {
'solution': {p: x[p].varValue for p in products},
'VaR': var.varValue,
'CVaR': value(model.objective)
}
Multi-Stage Stochastic Programming
Scenario Tree
Stage 1 → Stage 2 → Stage 3
x₁ → x₂(ξ₁) → x₃(ξ₁,ξ₂)
→ x₂(ξ₂) → x₃(ξ₂,ξ₃)
Dynamic Programming Approach
def multistage_inventory_dp(T, scenarios_per_stage):
"""
Multi-stage inventory control with dynamic programming
T: number of stages
scenarios_per_stage: number of scenarios at each stage
"""
# Value functions
V = [{} for _ in range(T+1)]
V[T] = {state: 0 for state in states} # Terminal value
# Backward recursion
for t in range(T-1, -1, -1):
for state in states:
min_cost = float('inf')
best_action = None
for action in actions:
# Expected cost-to-go
expected_cost = 0
for scenario in scenarios[t]:
next_state = transition(state, action, scenario)
prob = scenario['probability']
immediate_cost = cost(state, action, scenario)
future_cost = V[t+1][next_state]
expected_cost += prob * (immediate_cost + future_cost)
if expected_cost < min_cost:
min_cost = expected_cost
best_action = action
V[t][state] = min_cost
return V
Tools & Libraries
Python:
scipy.stats: probability distributionsnumpy: random samplingpulp/pyomo: stochastic programming formulationSALib: sensitivity analysis
Specialized:
PySP (Pyomo): stochastic programming extensionStochOptim.jl(Julia): stochastic optimization
Commercial:
CPLEX Stochastic SolverGurobi Multi-Scenario
Related Skills
- optimization-modeling: deterministic optimization
- demand-forecasting: uncertainty modeling
- inventory-optimization: stochastic inventory
- risk-mitigation: risk management
- scenario-planning: scenario generation
Related skills
More from kishorkukreja/awesome-supply-chain
procurement-optimization
When the user wants to optimize procurement decisions, allocate orders across suppliers, or determine optimal order quantities. Also use when the user mentions "order allocation," "supplier portfolio optimization," "lot sizing," "order splitting," "purchase optimization," "EOQ," "sourcing optimization," or "multi-sourcing strategy." For supplier selection, see supplier-selection. For spend analysis, see spend-analysis.
77replenishment-strategy
When the user wants to design or optimize replenishment strategies, determine replenishment policies, or improve inventory flow between locations. Also use when the user mentions "inventory replenishment," "stock replenishment," "min-max inventory," "DRP," "auto-replenishment," "vendor-managed inventory," "forward pick replenishment," or "retail store replenishment." For safety stock calculations, see inventory-optimization. For multi-echelon networks, see multi-echelon-inventory.
37inventory-optimization
When the user wants to optimize inventory levels, calculate safety stock, determine reorder points, or minimize inventory costs. Also use when the user mentions "inventory management," "safety stock," "EOQ," "reorder point," "service level," "stockout prevention," "ABC analysis," "inventory turns," or "working capital reduction." For warehouse slotting, see warehouse-slotting-optimization. For multi-echelon systems, see multi-echelon-inventory.
34supplier-selection
When the user wants to evaluate suppliers, select vendors, or perform supplier scoring and qualification. Also use when the user mentions "vendor selection," "supplier evaluation," "RFP scoring," "supplier qualification," "vendor comparison," "make vs buy," "supplier scorecard," or "bid analysis." For ongoing supplier risk monitoring, see supplier-risk-management. For contract negotiation, see contract-management.
32pharmaceutical-supply-chain
When the user wants to optimize pharmaceutical supply chains, manage cold chain logistics, ensure regulatory compliance, or implement serialization. Also use when the user mentions "pharma supply chain," "GMP compliance," "cold chain," "drug serialization," "clinical trials logistics," "pharmaceutical distribution," "good distribution practices," "GDP," "drug safety," or "pharmaceutical quality." For general healthcare, see hospital-logistics. For clinical trials specifically, see clinical-trial-logistics.
30retail-replenishment
When the user wants to optimize retail store replenishment, calculate reorder points for stores, or manage continuous replenishment. Also use when the user mentions "store replenishment," "auto-replenishment," "min-max inventory," "store orders," "DC to store," "continuous replenishment," or "vendor-managed inventory (VMI)." For initial allocation, see retail-allocation. For DC operations, see warehouse-design.
29