stochastic-optimization

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:

1. Uncertainty Characteristics

   • What parameters are uncertain? (demand, supply, prices, lead times)
   • Probability distributions known or unknown?
   • Historical data available?
   • Uncertainty independent or correlated?

2. Decision Structure

   • Single-stage or multi-stage decisions?
   • Which decisions made before/after uncertainty reveals?
   • Recourse actions available?
   • Decision frequency?

3. Risk Attitude

   • Risk-neutral (expected value) or risk-averse?
   • Preferred risk measure? (CVaR, variance, worst-case)
   • Service level requirements?
   • Budget/capacity constraints?

4. 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()