kishorkukreja/awesome-supply-chain/newsvendor-problem

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

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

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)