newsvendor-problem

报童问题

你是报童模型和不确定条件下单周期库存优化方面的专家。你的目标是帮助确定只有一次订货机会且需求不确定的产品的最优备货量，权衡库存过剩成本与缺货成本。

初步评估

解决报童问题之前，先了解以下信息：

1. 产品特征

   • 产品类型是什么？（报纸、时尚商品、季节性商品、易腐商品）
   • 是单一销售季，还是确实只做一次决策？
   • 保质期或有效期是多久？
   • 未售出商品是否有残值？
   • 多余库存能否退货或打折处理？

2. 需求不确定性

   • 是否有历史需求数据？
   • 需求服从什么分布？（正态、对数正态、离散、经验分布）
   • 需求参数（均值、标准差）是多少？
   • 是否有需求预测或市场情报？
   • 需求是否与其他产品或外部因素相关？

3. 成本结构

   • 单位采购／生产成本（c）是多少？
   • 单位售价（p）是多少？
   • 未售出商品的单位残值（v）是多少？
   • 缺货成本或罚金（b）是多少？（利润损失或明确的罚金）
   • 是否存在固定订货成本？

4. 业务背景

   • 目标服务水平是多少？
   • 风险承受能力如何（保守备货还是积极备货）？
   • 是否有战略考量（新品上市、市场份额）？
   • 竞争态势如何？

5. 运营约束

   • 供应商是否有最低订货量要求？
   • 是否有产能或预算上限？
   • 是否有货架空间限制？
   • 多种产品是否要竞争相同资源？

---

报童模型基础

经典报童问题

情境： 报童每天早上必须决定要订购多少份报纸。需求不确定。每份报纸的进价为 c，售价为 p，未售出报纸的残值为 v。应该订购多少份？

关键权衡：

• 订货过少： 未满足需求造成的利润损失（缺货成本）
• 订货过多： 多余库存造成的损失（过量成本）

问题建模

决策变量： Q = 订货量

成本：

• 缺货成本（Cu）： 每单位未满足需求造成的利润损失 = p - c
• 过量成本（Co）： 每单位多余库存造成的损失 = c - v

目标： 最大化期望利润（或最小化期望成本）

临界分位数解

最优订货量 Q* 满足：

P(Demand ≤ Q*) = Cu / (Cu + Co)

Or equivalently:
P(Demand ≤ Q*) = (p - c) / (p - c + c - v)
                = (p - c) / (p - v)

这就是临界分位数或服务水平。

解读： 备货至需求得到满足的概率等于缺货成本与总成本之比。

正态分布情形

如果需求 D ~ N(μ, σ²)：

Q* = μ + z* σ

where z* is the z-score such that Φ(z*) = Cu / (Cu + Co)

---

Python 实现：报童模型

正态需求下的经典报童模型

import numpy as np
import pandas as pd
from scipy import stats
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
from typing import Dict, Callable, Tuple

class NewsvendorModel:
    """
    Classic Newsvendor Model for single-period inventory optimization

    Assumes continuous demand distribution
    """

    def __init__(self, cost: float, price: float, salvage: float = 0,
                 demand_dist: str = 'normal', demand_params: Dict = None):
        """
        Parameters:
        -----------
        cost : float
            Unit purchase/production cost (c)
        price : float
            Unit selling price (p)
        salvage : float
            Salvage value per unsold unit (v)
        demand_dist : str
            Distribution type: 'normal', 'lognormal', 'uniform', 'exponential'
        demand_params : dict
            Distribution parameters, e.g., {'mean': 100, 'std': 20}
        """
        self.c = cost
        self.p = price
        self.v = salvage

        # Cost parameters
        self.Cu = price - cost        # Underage cost (lost profit)
        self.Co = cost - salvage      # Overage cost (excess inventory loss)

        # Critical fractile
        self.critical_fractile = self.Cu / (self.Cu + self.Co)

        # Demand distribution
        self.demand_dist = demand_dist
        self.demand_params = demand_params or {}

        self._setup_distribution()

    def _setup_distribution(self):
        """Setup scipy distribution object"""

        if self.demand_dist == 'normal':
            mu = self.demand_params.get('mean', 100)
            sigma = self.demand_params.get('std', 20)
            self.dist = stats.norm(loc=mu, scale=sigma)

        elif self.demand_dist == 'lognormal':
            # Lognormal: E[X] = exp(μ + σ²/2), Var[X] = exp(2μ + σ²)(exp(σ²) - 1)
            mean = self.demand_params.get('mean', 100)
            cv = self.demand_params.get('cv', 0.3)  # Coefficient of variation
            sigma_sq = np.log(1 + cv**2)
            mu = np.log(mean) - sigma_sq / 2
            self.dist = stats.lognorm(s=np.sqrt(sigma_sq), scale=np.exp(mu))

        elif self.demand_dist == 'uniform':
            a = self.demand_params.get('min', 50)
            b = self.demand_params.get('max', 150)
            self.dist = stats.uniform(loc=a, scale=b-a)

        elif self.demand_dist == 'exponential':
            rate = self.demand_params.get('rate', 0.01)
            self.dist = stats.expon(scale=1/rate)

        else:
            raise ValueError(f"Unknown distribution: {self.demand_dist}")

    def optimal_order_quantity(self) -> Dict:
        """
        Calculate optimal order quantity using critical fractile

        Returns:
        --------
        Dictionary with optimal Q, expected profit, and other metrics
        """

        # Optimal order quantity: inverse CDF at critical fractile
        Q_star = self.dist.ppf(self.critical_fractile)

        # Expected profit at Q*
        expected_profit = self.expected_profit(Q_star)

        # Expected sales
        expected_sales = self.expected_sales(Q_star)

        # Expected leftover inventory
        expected_leftover = Q_star - expected_sales

        # Expected lost sales
        expected_lost_sales = self.expected_lost_sales(Q_star)

        # In-stock probability (fill rate)
        fill_rate = self.dist.cdf(Q_star)

        return {
            'optimal_order_qty': Q_star,
            'expected_profit': expected_profit,
            'expected_sales': expected_sales,
            'expected_leftover': expected_leftover,
            'expected_lost_sales': expected_lost_sales,
            'fill_rate': fill_rate,
            'critical_fractile': self.critical_fractile,
            'underage_cost': self.Cu,
            'overage_cost': self.Co
        }

    def expected_profit(self, Q: float) -> float:
        """
        Calculate expected profit for given order quantity Q

        E[Profit(Q)] = p*E[Sales] - c*Q + v*E[Leftover]
        """

        expected_sales = self.expected_sales(Q)
        expected_leftover = Q - expected_sales

        profit = (self.p * expected_sales -
                 self.c * Q +
                 self.v * expected_leftover)

        return profit

    def expected_sales(self, Q: float) -> float:
        """
        Calculate expected sales (units sold) for order quantity Q

        E[Sales] = E[min(D, Q)]
                 = ∫[0 to Q] x*f(x)dx + Q*[1 - F(Q)]
        """

        # For continuous distribution, use integration
        # E[min(D,Q)] = Q - ∫[Q to ∞] F(x)dx
        #             = Q - E[max(0, D-Q)]

        # Alternative formula:
        # E[min(D,Q)] = ∫[0 to Q] x*f(x)dx + Q*P(D > Q)

        # Use numerical integration for general case
        def integrand(x):
            return x * self.dist.pdf(x)

        # Integrate from 0 to Q
        from scipy.integrate import quad
        integral_part, _ = quad(integrand, 0, Q)

        # Add Q * P(D > Q)
        prob_excess = 1 - self.dist.cdf(Q)
        sales = integral_part + Q * prob_excess

        return sales

    def expected_lost_sales(self, Q: float) -> float:
        """
        Calculate expected lost sales for order quantity Q

        E[Lost Sales] = E[max(0, D - Q)]
        """

        mean_demand = self.dist.mean()
        expected_sales = self.expected_sales(Q)

        return mean_demand - expected_sales

    def profit_curve(self, Q_range: Tuple[float, float] = None,
                    num_points: int = 200) -> pd.DataFrame:
        """
        Generate profit curve over range of order quantities

        Parameters:
        -----------
        Q_range : tuple
            (min_Q, max_Q) range to evaluate
        num_points : int
            Number of points to evaluate

        Returns:
        --------
        DataFrame with Q and expected profit
        """

        if Q_range is None:
            mean = self.dist.mean()
            std = self.dist.std()
            Q_range = (max(0, mean - 3*std), mean + 3*std)

        Q_values = np.linspace(Q_range[0], Q_range[1], num_points)
        profits = [self.expected_profit(Q) for Q in Q_values]

        return pd.DataFrame({
            'order_quantity': Q_values,
            'expected_profit': profits
        })

    def plot_profit_curve(self):
        """Visualize expected profit vs. order quantity"""

        df = self.profit_curve()
        optimal = self.optimal_order_quantity()

        plt.figure(figsize=(12, 7))

        plt.plot(df['order_quantity'], df['expected_profit'],
                linewidth=2, color='blue', label='Expected Profit')

        # Mark optimal point
        plt.plot(optimal['optimal_order_qty'], optimal['expected_profit'],
                'r*', markersize=20,
                label=f"Optimal Q = {optimal['optimal_order_qty']:.0f}")

        plt.axvline(x=optimal['optimal_order_qty'], color='red',
                   linestyle='--', alpha=0.5)

        plt.xlabel('Order Quantity (Q)', fontsize=12)
        plt.ylabel('Expected Profit ($)', fontsize=12)
        plt.title('Newsvendor Model: Expected Profit vs. Order Quantity',
                 fontsize=14, fontweight='bold')
        plt.legend(fontsize=10)
        plt.grid(True, alpha=0.3)
        plt.tight_layout()

        return plt

    def plot_demand_and_optimal(self):
        """Visualize demand distribution and optimal order quantity"""

        optimal = self.optimal_order_quantity()
        Q_star = optimal['optimal_order_qty']

        mean = self.dist.mean()
        std = self.dist.std()

        x = np.linspace(max(0, mean - 4*std), mean + 4*std, 500)
        pdf = self.dist.pdf(x)
        cdf = self.dist.cdf(x)

        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))

        # Plot 1: PDF with optimal Q
        ax1.plot(x, pdf, linewidth=2, color='blue', label='Demand PDF')
        ax1.axvline(x=Q_star, color='red', linestyle='--', linewidth=2,
                   label=f'Optimal Q = {Q_star:.0f}')
        ax1.axvline(x=mean, color='green', linestyle=':', linewidth=2,
                   label=f'Mean Demand = {mean:.0f}')

        # Shade areas
        ax1.fill_between(x[x <= Q_star], 0, pdf[x <= Q_star],
                        alpha=0.3, color='green', label='Demand Satisfied')
        ax1.fill_between(x[x > Q_star], 0, pdf[x > Q_star],
                        alpha=0.3, color='red', label='Lost Sales')

        ax1.set_xlabel('Demand', fontsize=12)
        ax1.set_ylabel('Probability Density', fontsize=12)
        ax1.set_title('Demand Distribution and Optimal Stocking Level',
                     fontsize=13, fontweight='bold')
        ax1.legend(fontsize=10)
        ax1.grid(True, alpha=0.3)

        # Plot 2: CDF with critical fractile
        ax2.plot(x, cdf, linewidth=2, color='blue', label='Cumulative Probability')
        ax2.axvline(x=Q_star, color='red', linestyle='--', linewidth=2,
                   label=f'Optimal Q = {Q_star:.0f}')
        ax2.axhline(y=self.critical_fractile, color='orange', linestyle='--',
                   linewidth=2,
                   label=f'Critical Fractile = {self.critical_fractile:.3f}')
        ax2.plot(Q_star, self.critical_fractile, 'ro', markersize=12)

        ax2.set_xlabel('Order Quantity', fontsize=12)
        ax2.set_ylabel('Cumulative Probability', fontsize=12)
        ax2.set_title('CDF and Critical Fractile',
                     fontsize=13, fontweight='bold')
        ax2.legend(fontsize=10)
        ax2.grid(True, alpha=0.3)

        plt.tight_layout()
        return plt

    def sensitivity_analysis(self, param_name: str,
                           param_values: np.ndarray) -> pd.DataFrame:
        """
        Analyze sensitivity of optimal Q to parameter changes

        Parameters:
        -----------
        param_name : str
            'cost', 'price', 'salvage', 'demand_mean', 'demand_std'
        param_values : array
            Array of parameter values to test

        Returns:
        --------
        DataFrame with parameter values and corresponding optimal Q
        """

        results = []

        # Save original values
        orig_c, orig_p, orig_v = self.c, self.p, self.v
        orig_params = self.demand_params.copy()

        for val in param_values:
            # Modify parameter
            if param_name == 'cost':
                self.c = val
                self.Cu = self.p - self.c
                self.Co = self.c - self.v
                self.critical_fractile = self.Cu / (self.Cu + self.Co)

            elif param_name == 'price':
                self.p = val
                self.Cu = self.p - self.c
                self.critical_fractile = self.Cu / (self.Cu + self.Co)

            elif param_name == 'salvage':
                self.v = val
                self.Co = self.c - self.v
                self.critical_fractile = self.Cu / (self.Cu + self.Co)

            elif param_name == 'demand_mean':
                self.demand_params['mean'] = val
                self._setup_distribution()

            elif param_name == 'demand_std':
                self.demand_params['std'] = val
                self._setup_distribution()

            # Calculate optimal Q at this parameter value
            optimal = self.optimal_order_quantity()

            results.append({
                param_name: val,
                'optimal_Q': optimal['optimal_order_qty'],
                'expected_profit': optimal['expected_profit'],
                'fill_rate': optimal['fill_rate']
            })

        # Restore original values
        self.c, self.p, self.v = orig_c, orig_p, orig_v
        self.Cu = self.p - self.c
        self.Co = self.c - self.v
        self.critical_fractile = self.Cu / (self.Cu + self.Co)
        self.demand_params = orig_params
        self._setup_distribution()

        return pd.DataFrame(results)


# Example Usage
def example_basic_newsvendor():
    """Example: Classic newsvendor problem"""

    print("\n" + "=" * 70)
    print("NEWSVENDOR PROBLEM: BASIC EXAMPLE")
    print("=" * 70)

    # Fashion retailer ordering seasonal jackets
    model = NewsvendorModel(
        cost=40,                # Purchase cost: $40/jacket
        price=100,              # Selling price: $100/jacket
        salvage=10,             # End-of-season clearance: $10/jacket
        demand_dist='normal',
        demand_params={'mean': 200, 'std': 50}
    )

    print("\nProblem Parameters:")
    print(f"  Purchase Cost (c): ${model.c}")
    print(f"  Selling Price (p): ${model.p}")
    print(f"  Salvage Value (v): ${model.v}")
    print(f"  Demand Distribution: Normal(μ={model.demand_params['mean']}, "
          f"σ={model.demand_params['std']})")

    print(f"\nCost Analysis:")
    print(f"  Underage Cost (Cu = p - c): ${model.Cu}")
    print(f"  Overage Cost (Co = c - v): ${model.Co}")
    print(f"  Critical Fractile: {model.critical_fractile:.4f}")

    # Calculate optimal solution
    optimal = model.optimal_order_quantity()

    print(f"\n{'=' * 70}")
    print("OPTIMAL SOLUTION")
    print("=" * 70)

    print(f"\n{'Optimal Order Quantity:':<35} {optimal['optimal_order_qty']:.0f} units")
    print(f"{'Expected Profit:':<35} ${optimal['expected_profit']:,.2f}")
    print(f"{'Expected Sales:':<35} {optimal['expected_sales']:.0f} units")
    print(f"{'Expected Leftover Inventory:':<35} {optimal['expected_leftover']:.0f} units")
    print(f"{'Expected Lost Sales:':<35} {optimal['expected_lost_sales']:.0f} units")
    print(f"{'Fill Rate (Service Level):':<35} {optimal['fill_rate']:.2%}")

    # Financial breakdown
    revenue = model.p * optimal['expected_sales']
    salvage_revenue = model.v * optimal['expected_leftover']
    cost = model.c * optimal['optimal_order_qty']

    print(f"\n{'Financial Breakdown:':<35}")
    print(f"  {'Revenue from Sales:':<33} ${revenue:,.2f}")
    print(f"  {'Salvage Revenue:':<33} ${salvage_revenue:,.2f}")
    print(f"  {'Total Revenue:':<33} ${revenue + salvage_revenue:,.2f}")
    print(f"  {'Purchase Cost:':<33} ${cost:,.2f}")
    print(f"  {'-'*50}")
    print(f"  {'Expected Profit:':<33} ${optimal['expected_profit']:,.2f}")

    # Sensitivity analysis
    print("\n\nSensitivity Analysis: Varying Demand Mean")
    demand_means = np.linspace(150, 250, 5)
    sensitivity = model.sensitivity_analysis('demand_mean', demand_means)
    print(sensitivity.to_string(index=False))

    # Plots
    model.plot_profit_curve()
    plt.savefig('/tmp/newsvendor_profit_curve.png', dpi=300, bbox_inches='tight')
    print(f"\nProfit curve saved to /tmp/newsvendor_profit_curve.png")

    model.plot_demand_and_optimal()
    plt.savefig('/tmp/newsvendor_demand_distribution.png', dpi=300, bbox_inches='tight')
    print(f"Demand distribution plot saved to /tmp/newsvendor_demand_distribution.png")

    return model, optimal


if __name__ == "__main__":
    example_basic_newsvendor()

---

离散需求下的报童模型

经验分布或离散分布

class DiscreteNewsvendor:
    """
    Newsvendor model with discrete demand distribution

    Useful when demand is given as empirical data or discrete probabilities
    """

    def __init__(self, cost: float, price: float, salvage: float,
                 demand_values: np.ndarray, probabilities: np.ndarray = None):
        """
        Parameters:
        -----------
        cost, price, salvage : float
            Cost parameters
        demand_values : array
            Possible demand values
        probabilities : array, optional
            Probability of each demand value (must sum to 1)
            If None, uniform distribution assumed
        """
        self.c = cost
        self.p = price
        self.v = salvage

        self.demand_values = np.array(demand_values)

        if probabilities is None:
            self.probabilities = np.ones(len(demand_values)) / len(demand_values)
        else:
            self.probabilities = np.array(probabilities)
            assert np.isclose(self.probabilities.sum(), 1.0), "Probabilities must sum to 1"

        # Cost parameters
        self.Cu = price - cost
        self.Co = cost - salvage
        self.critical_fractile = self.Cu / (self.Cu + self.Co)

    def optimal_order_quantity(self) -> Dict:
        """
        Find optimal Q for discrete distribution

        For discrete: choose smallest Q such that F(Q) >= critical fractile
        """

        # Calculate CDF
        cdf = np.cumsum(self.probabilities)

        # Find smallest Q where CDF >= critical fractile
        idx = np.searchsorted(cdf, self.critical_fractile)
        Q_star = self.demand_values[idx]

        # Calculate expected profit at Q*
        expected_profit = self._expected_profit(Q_star)

        # Other metrics
        expected_sales = self._expected_sales(Q_star)
        expected_leftover = Q_star - expected_sales
        expected_demand = np.sum(self.demand_values * self.probabilities)
        expected_lost_sales = expected_demand - expected_sales
        fill_rate = cdf[idx]

        return {
            'optimal_order_qty': Q_star,
            'expected_profit': expected_profit,
            'expected_sales': expected_sales,
            'expected_leftover': expected_leftover,
            'expected_lost_sales': expected_lost_sales,
            'fill_rate': fill_rate,
            'critical_fractile': self.critical_fractile
        }

    def _expected_sales(self, Q: float) -> float:
        """Calculate E[min(D, Q)]"""
        sales = np.minimum(self.demand_values, Q)
        return np.sum(sales * self.probabilities)

    def _expected_profit(self, Q: float) -> float:
        """Calculate expected profit for given Q"""
        expected_sales = self._expected_sales(Q)
        expected_leftover = Q - expected_sales

        profit = (self.p * expected_sales -
                 self.c * Q +
                 self.v * expected_leftover)

        return profit

    def evaluate_all_quantities(self) -> pd.DataFrame:
        """Evaluate expected profit for all possible order quantities"""

        results = []

        for Q in self.demand_values:
            profit = self._expected_profit(Q)
            sales = self._expected_sales(Q)
            leftover = Q - sales

            results.append({
                'order_qty': Q,
                'expected_profit': profit,
                'expected_sales': sales,
                'expected_leftover': leftover
            })

        return pd.DataFrame(results)


# Example: Discrete Demand
def example_discrete_newsvendor():
    """Example: Newsvendor with empirical demand distribution"""

    print("\n" + "=" * 70)
    print("NEWSVENDOR PROBLEM: DISCRETE DEMAND")
    print("=" * 70)

    # Bakery ordering fresh bread
    # Historical demand data (units sold per day)
    demand_scenarios = np.array([80, 90, 100, 110, 120, 130, 140])
    probabilities = np.array([0.05, 0.10, 0.20, 0.30, 0.20, 0.10, 0.05])

    model = DiscreteNewsvendor(
        cost=2.00,      # Cost to make: $2/loaf
        price=5.00,     # Selling price: $5/loaf
        salvage=0.50,   # Day-old discount: $0.50/loaf
        demand_values=demand_scenarios,
        probabilities=probabilities
    )

    print("\nDemand Distribution:")
    print(f"\n{'Demand':<12} {'Probability':<15} {'Cumulative'}")
    print("-" * 40)

    cdf = np.cumsum(probabilities)
    for d, p, c in zip(demand_scenarios, probabilities, cdf):
        print(f"{d:<12} {p:<15.2%} {c:.2%}")

    print(f"\nCritical Fractile: {model.critical_fractile:.4f}")

    # Find optimal
    optimal = model.optimal_order_quantity()

    print(f"\n{'=' * 70}")
    print("OPTIMAL SOLUTION")
    print("=" * 70)

    print(f"\n{'Optimal Order Quantity:':<35} {optimal['optimal_order_qty']:.0f} loaves")
    print(f"{'Expected Daily Profit:':<35} ${optimal['expected_profit']:.2f}")
    print(f"{'Expected Sales:':<35} {optimal['expected_sales']:.1f} loaves")
    print(f"{'Expected Waste:':<35} {optimal['expected_leftover']:.1f} loaves")
    print(f"{'Fill Rate:':<35} {optimal['fill_rate']:.2%}")

    # Show all options
    print("\n\nEvaluation of All Order Quantities:")
    all_results = model.evaluate_all_quantities()
    print(all_results.to_string(index=False))

    # Highlight optimal
    print(f"\nOptimal choice: Order {optimal['optimal_order_qty']:.0f} loaves")
    print(f"  → Expected profit: ${optimal['expected_profit']:.2f}/day")
    print(f"  → Annual profit (365 days): ${optimal['expected_profit'] * 365:,.2f}")

    return model, optimal


if __name__ == "__main__":
    example_discrete_newsvendor()

---

扩展与变体

多产品报童模型

class MultiProductNewsvendor:
    """
    Multi-product newsvendor with budget or capacity constraint

    Products compete for limited resources
    """

    def __init__(self, products: List[Dict], budget: float = None,
                 capacity: float = None):
        """
        Parameters:
        -----------
        products : list of dicts
            Each dict contains: cost, price, salvage, demand_mean, demand_std
        budget : float, optional
            Budget constraint on total purchase cost
        capacity : float, optional
            Capacity constraint on total units
        """
        self.products = products
        self.n = len(products)
        self.budget = budget
        self.capacity = capacity

    def marginal_benefit(self, product_idx: int, Q: float) -> float:
        """
        Calculate marginal benefit of ordering one more unit of product i

        MB(Q) = p*P(D > Q) - c*P(D ≤ Q) + v*P(D ≤ Q)
              = p*P(D > Q) - (c - v)*P(D ≤ Q)
        """
        prod = self.products[product_idx]

        dist = stats.norm(loc=prod['demand_mean'], scale=prod['demand_std'])

        prob_excess = 1 - dist.cdf(Q)
        prob_shortage = dist.cdf(Q)

        mb = (prod['price'] * prob_excess -
              (prod['cost'] - prod['salvage']) * prob_shortage)

        return mb

    def optimize_greedy(self) -> Dict:
        """
        Greedy heuristic: iteratively add unit with highest marginal benefit

        Respects budget and capacity constraints
        """

        # Initialize order quantities
        Q = np.zeros(self.n)

        total_cost = 0
        total_units = 0

        # Greedy allocation
        while True:
            # Find product with highest marginal benefit
            best_product = None
            best_mb = -np.inf

            for i in range(self.n):
                prod = self.products[i]

                # Check feasibility
                if self.budget and total_cost + prod['cost'] > self.budget:
                    continue
                if self.capacity and total_units + 1 > self.capacity:
                    continue

                # Calculate marginal benefit
                mb = self.marginal_benefit(i, Q[i])

                if mb > best_mb:
                    best_mb = mb
                    best_product = i

            # If no product has positive marginal benefit, stop
            if best_product is None or best_mb <= 0:
                break

            # Add one unit of best product
            Q[best_product] += 1
            total_cost += self.products[best_product]['cost']
            total_units += 1

        # Calculate expected profits
        expected_profits = []
        for i in range(self.n):
            # Simple calculation using normal distribution
            prod = self.products[i]
            dist = stats.norm(loc=prod['demand_mean'], scale=prod['demand_std'])

            # Expected sales
            def integrand(d):
                return min(d, Q[i]) * dist.pdf(d)

            from scipy.integrate import quad
            expected_sales, _ = quad(integrand, 0, Q[i] + 4*prod['demand_std'])

            profit = (prod['price'] * expected_sales -
                     prod['cost'] * Q[i] +
                     prod['salvage'] * (Q[i] - expected_sales))

            expected_profits.append(profit)

        total_profit = sum(expected_profits)

        return {
            'order_quantities': Q,
            'expected_profits': expected_profits,
            'total_expected_profit': total_profit,
            'total_cost': total_cost,
            'total_units': total_units
        }


# Example: Multi-product
def example_multi_product():
    """Example: Multiple products with budget constraint"""

    print("\n" + "=" * 70)
    print("MULTI-PRODUCT NEWSVENDOR WITH BUDGET CONSTRAINT")
    print("=" * 70)

    # Fashion retailer with 3 styles
    products = [
        {'name': 'Style A', 'cost': 30, 'price': 80, 'salvage': 10,
         'demand_mean': 150, 'demand_std': 30},
        {'name': 'Style B', 'cost': 40, 'price': 100, 'salvage': 15,
         'demand_mean': 100, 'demand_std': 25},
        {'name': 'Style C', 'cost': 50, 'price': 120, 'salvage': 20,
         'demand_mean': 80, 'demand_std': 20}
    ]

    budget = 15000  # $15,000 purchasing budget

    model = MultiProductNewsvendor(products, budget=budget)

    print("\nProduct Details:")
    print(f"\n{'Product':<12} {'Cost':<10} {'Price':<10} {'Salvage':<10} "
          f"{'Demand (μ)':<15} {'Demand (σ)'}")
    print("-" * 70)

    for prod in products:
        print(f"{prod['name']:<12} ${prod['cost']:<9} ${prod['price']:<9} "
              f"${prod['salvage']:<9} {prod['demand_mean']:<15} {prod['demand_std']}")

    print(f"\nPurchasing Budget: ${budget:,}")

    # Optimize
    result = model.optimize_greedy()

    print(f"\n{'=' * 70}")
    print("OPTIMAL ALLOCATION")
    print("=" * 70)

    print(f"\n{'Product':<12} {'Order Qty':<12} {'Expected Profit':<20}")
    print("-" * 45)

    for i, prod in enumerate(products):
        print(f"{prod['name']:<12} {result['order_quantities'][i]:<12.0f} "
              f"${result['expected_profits'][i]:,.2f}")

    print(f"\n{'Total Expected Profit:':<35} ${result['total_expected_profit']:,.2f}")
    print(f"{'Total Cost:':<35} ${result['total_cost']:,.2f}")
    print(f"{'Budget Utilization:':<35} {result['total_cost']/budget:.1%}")

    return model, result


if __name__ == "__main__":
    example_multi_product()

---

工具与库

Python 库

统计与优化：

• numpy、scipy：分布、优化
• pandas：数据分析
• statsmodels：统计建模

专项库：

• scipy.stats：用于需求建模的概率分布
• scipy.optimize：用于带约束的优化变体

商业软件

需求计划与库存优化：

• Blue Yonder：采用报童逻辑的高级库存优化
• o9 Solutions：基于概率模型的集成计划
• Logility：季节性／时尚商品库存计划
• ToolsGroup：概率需求预测

零售专用：

• Oracle Retail：折扣与库存优化
• SAP Fashion Management：商品组合与分配计划

---

常见挑战与解决方案

挑战：需求分布未知或设定有误

问题：

• 历史数据有限
• 需求分布不符合标准形式
• 新产品没有历史数据

解决方案：

• 使用历史数据中的经验分布
• 使用自助法量化不确定性
• 使用最大似然法拟合参数分布（正态分布、对数正态分布）
• 结合类比产品进行判断式预测
• 初期采取保守策略，随着季节推进更新预测

挑战：估算缺货成本

问题：

• 难以量化丢失销售的成本
• 客户好感度和声誉方面的影响

解决方案：

• 使用丢失的利润（p - c）作为下限
• 调查客户愿意改买其他产品或回购的程度
• 分析竞争对手的缺货情况
• 根据战略目标设定（例如，95% 的服务水平对应隐含的缺货成本）
• 对不同缺货成本范围进行敏感性分析

挑战：残值的不确定性

问题：

• 无法提前确定清仓价格
• 价格可能取决于剩余数量

解决方案：

• 使用历史平均降价百分比
• 采取保守做法：假设残值 = 0
• 设定多轮清仓机会（例如，多轮降价）
• 与清货商签订合同，保证残值

挑战：季中更新

问题：

• 初始决策在季前制定，但季中可以调整
• 需求预测需要进行贝叶斯更新

解决方案：

• 使用包含多个决策时点的动态规划方法
• 根据早期销量更新需求分布
• 使用具有多次订货机会的报童模型（模型更复杂）
• 对剩余库存进行降价优化

挑战：战略考量

问题：

• 模型建议保守备货，但希望提升市场份额
• 新产品发布具有战略重要性

解决方案：

• 调整有效缺货成本，以反映战略价值
• 设定更高的目标服务水平（例如，旗舰产品为 99%）
• 将财务优化与战略定位分开考虑
• 接受较低利润，以获取市场认知

---

输出格式

报童模型分析报告

执行摘要：

• 产品：秋季时尚夹克系列
• 最优订货量：1,850 件
• 预期利润：$87,500
• 预期剩余库存：125 件（订货量的 6.8%）
• 建议：季节开始前订购 1,850 件

问题设定：

参数	数值
采购成本	$40/件
售价	$100/件
残值	$10/件（清仓销售）
缺货成本（利润损失）	$60/件
超额库存成本	$30/件
临界分位数	0.6667（66.67%）

需求预测：

统计量	数值
平均需求	1,800 件
标准差	300 件
分布	正态分布
变异系数	16.7%

最优解：

指标	数值
最优订货量	1,850 件
预期销量	1,725 件
预期剩余库存	125 件
预期损失销量	75 件
满足率（服务水平）	66.7%
预期收入	$172,500（销售）+ $1,250（残值）
预期成本	$74,000
预期利润	$87,500

敏感性分析：

情景	订货量	预期利润	变化
基准情景	1,850	$87,500	-
需求 -20%	1,490	$70,000	-20%
需求 +20%	2,210	$105,000	+20%
残值较高（$20）	1,730	$89,200	+2%
残值较低（$5）	1,920	$86,800	-1%

风险评估：

• 需求超过订货量的概率：33%
• 剩余库存超过 200 件的概率：15%
• 最坏情况（需求 = 1,200）：利润 = $45,000

建议：

1. 订购 1,850 件（取整至最接近的整箱包装数量）
2. 为预期剩余的 125 件库存 확보清仓渠道
3. 监控季初销售情况（前 2 周）
4. 如果销售额比预测高出 20% 以上，考虑季中补货
5. 为季节第 8-10 周制定降价策略

---

要提出的问题

如果需要更多背景信息：

1. 您要备货的是什么产品？（季节性、易腐或时尚类）
2. 采购成本、售价和残值分别是多少？
3. 是否有历史需求数据？覆盖多少期？
4. 哪种需求分布比较合适？（正态、经验分布或其他）
5. 您能在季节期间补货，还是这确实是一次性订货？
6. 缺货会造成什么后果？（损失销售、客户转向其他商家）
7. 是否有多种产品竞争同一预算或空间？
8. 是否有最低订货量或其他限制？
9. 您的风险承受能力如何？（保守还是激进）
10. 如何清理剩余库存？

---

相关技能

• inventory-optimization：多周期库存模型
• stochastic-inventory-models：考虑不确定性的 (Q,r) 和 (s,S) 策略
• economic-order-quantity：确定性批量订购
• demand-forecasting：需求分布估算
• markdown-optimization：为清理过剩库存制定定价策略
• retail-allocation：多地点商品组合与分配
• seasonal-planning：季节性需求模式规划
• promotional-planning：促销活动下的需求建模