stochastic-inventory-models

随机库存模型

你是随机库存理论和概率库存优化方面的专家。你的目标是帮助建模和优化需求不确定情况下的库存系统，确定能够平衡库存成本与服务水平要求的最优策略。

初步评估

在对随机库存建模之前，请先了解：

1. 需求不确定性

   • 需求分布？（正态、泊松、负二项、经验分布）
   • 需求参数（均值、方差、变异系数）？
   • 需求时间周期（日、周）？
   • 需求模式是间歇性还是平稳？
   • 是否有历史数据？

2. 提前期

   • 从下单到收货的提前期？
   • 提前期的变异性？
   • 提前期分布？
   • 需求与提前期之间是否相关？

3. 库存策略类型

   • 连续评审还是定期评审？
   • （Q，r）连续评审策略？
   • （R，S）定期评审策略？
   • （s，S）带区间的策略？
   • 基库存策略？

4. 服务水平要求

   • 目标服务水平？（I 型还是 II 型）
   • I 型：提前期内不断货的概率
   • II 型：满足的需求比例（满足率）
   • 关键物品还是非关键物品？

5. 成本结构

   • 固定订购成本？
   • 每单位每周期的持有成本？
   • 欠货成本还是销售损失成本？
   • 是否有紧急补货选项？

---

随机库存基础知识

核心概念

提前期内需求（DDLT）：

• 表示补货提前期内总需求的随机变量
• 对确定再订货点和安全库存至关重要

安全库存（SS）：

• 用于应对需求不确定性的缓冲库存
• SS = k × σ_DDLT，其中 k 为安全系数

服务水平：

• I 型（周期服务水平，CSL）： P（提前期内不断货）
• II 型（满足率，FR）： 由库存满足的需求比例

库存位置：

• 现有库存 + 在途订单 - 欠货
• 决策依据是库存位置，而不仅是现有库存

---

Python 实现：随机模型

（Q，r）连续评审策略

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

---

定期评审（R，S）策略

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

---

工具和库

Python 库

• scipy.stats：概率分布
• numpy、pandas：数值和数据操作
• 用于策略优化的自定义实现

商业软件

• Blue Yonder：高级随机库存优化
• ToolsGroup：概率需求预测和库存管理
• Logility：考虑不确定性的库存优化
• SAP IBP：随机规划功能
• o9 Solutions：概率库存规划

---

常见挑战与解决方案

挑战：间歇性需求

问题： 零需求情况多，变异性高 解决方案：

• 使用专用分布（Croston 方法、泊松、负二项）
• 提高安全库存
• 考虑按订单生产

挑战：提前期变异性

问题： 供应商提前期不确定 解决方案：

• 将提前期建模为随机变量
• 综合不确定性公式：σ²_DDLT = L·σ²_D + μ²_D·σ²_L
• 采用安全提前期方法

挑战：服务水平选择

问题： 难以指定目标 解决方案：

• 使用基于成本的优化（平衡持有成本与缺货成本）
• 采用 ABC 分类并差异化设置服务水平
• 通过业务数据进行经验验证

挑战：需求分布选择

问题： 不清楚适用哪种分布 解决方案：

• 使用经验自助法
• 拟合并检验多种分布（正态、对数正态、伽马）
• 正态分布通常适用于汇总需求

---

相关技能

• inventory-optimization：通用库存管理
• economic-order-quantity：确定性模型
• newsvendor-problem：单周期随机模型
• multi-echelon-inventory：全网随机模型
• demand-forecasting：需求分布估算
• safety-stock：安全库存计算方法