FC Inventory Optimizer

import numpy as np
from scipy import stats
from dataclasses import dataclass

@dataclass
class SKUParams:
    sku_id: str
    annual_demand: float          # D: forecasted annual units
    demand_std_weekly: float     # sigma: weekly demand std deviation
    lead_time_weeks: float        # L: vendor to FC transit time
    unit_cost: float               # C: landed cost per unit
    storage_cost_cuft_mo: float   # FC storage rate ($/cu ft/month)
    cube_per_unit: float          # cubic feet per unit
    supplier_moq: int             # minimum order quantity

def compute_holding_cost(p: SKUParams, capital_rate: float = 0.12,
                          obsolescence_rate: float = 0.03) -> float:
    capital = p.unit_cost * capital_rate
    storage = p.storage_cost_cuft_mo * p.cube_per_unit * 12
    obsolescence = p.unit_cost * obsolescence_rate
    return capital + storage + obsolescence

def eoq_with_moq(D: float, S: float, H: float, moq: int) -> int:
    raw_eoq = np.sqrt(2 * D * S / H)
    return max(moq, int(np.ceil(raw_eoq / moq) * moq))

def safety_stock(sigma_weekly: float, lead_time: float,
                   service_level: float = 0.95) -> int:
    z = stats.norm.ppf(service_level)  # 0.95 -> 1.645, 0.99 -> 2.326
    return int(np.ceil(z * sigma_weekly * np.sqrt(lead_time)))

def reorder_point(avg_weekly_demand: float, lead_time: float,
                    ss: int) -> int:
    return int(np.ceil(avg_weekly_demand * lead_time + ss))

def total_annual_cost(D: float, Q: int, S: float, H: float, ss: int) -> float:
    return (D / Q) * S + (Q / 2 + ss) * H

# Echo Dot 5th Gen at BOS5
params = SKUParams("B09B8V1LZ3", 52000, 85, 2, 22.99, 0.87, 0.18, 500)
H = compute_holding_cost(params)        # ~$5.65/unit/yr
S = 350.0                              # PO processing + freight + receiving
Q = eoq_with_moq(params.annual_demand, S, H, params.supplier_moq)
ss = safety_stock(params.demand_std_weekly, params.lead_time_weeks)
rop = reorder_point(params.annual_demand / 52, params.lead_time_weeks, ss)
tc = total_annual_cost(params.annual_demand, Q, S, H, ss)
print(f"EOQ={Q}, SS={ss}, ROP={rop}, TC=${tc:,.0f}")
# Output: EOQ=2000, SS=197, ROP=397, TC=$17,043

import numpy as np

class HoltWinters:
    """Triple exponential smoothing: level + trend + multiplicative seasonality."""
    def __init__(self, season_len: int = 52):
        self.m = season_len  # 52 weeks for weekly data

    def fit(self, y: np.ndarray, alpha=0.3, beta=0.05, gamma=0.15):
        m, n = self.m, len(y)
        level = np.mean(y[:m])
        trend = (np.mean(y[m:2*m]) - np.mean(y[:m])) / m
        seasonal = [y[i] - level for i in range(m)]
        self.fitted = []
        for t in range(n):
            forecast = (level + trend) * seasonal[t % m] if t >= m else y[t]
            self.fitted.append(forecast)
            prev = level
            level = alpha * (y[t] / seasonal[t % m]) + (1 - alpha) * (level + trend)
            trend = beta * (level - prev) + (1 - beta) * trend
            seasonal[t % m] = gamma * (y[t] / level) + (1 - gamma) * seasonal[t % m]
        self.level, self.trend, self.seasonal = level, trend, seasonal
        return self

    def forecast(self, horizon: int) -> np.ndarray:
        preds = []
        for h in range(1, horizon + 1):
            yhat = (self.level + h * self.trend) * self.seasonal[
                (len(self.fitted) + h) % self.m]
            preds.append(max(0, yhat))
        return np.array(preds)

    def mape(self, y: np.ndarray) -> float:
        fitted = np.array(self.fitted[self.m:])
        actual = y[self.m:]
        return float(np.mean(np.abs((actual - fitted) / actual)) * 100)

# Usage: fit on 104 weeks of history, forecast next 8 weeks for EOQ input
# demand_history = load_weekly_demand("B09B8V1LZ3", weeks=104)
# hw = HoltWinters(52).fit(demand_history, alpha=0.3, beta=0.05, gamma=0.15)
# next_8_weeks = hw.forecast(8)  # feed into EOQ annual_demand
# print(f"MAPE: {hw.mape(demand_history):.1f}%")  # target: <5%

import numpy as np
from scipy import stats
from dataclasses import dataclass
from typing import List

@dataclass
class Echelon:
    name: str
    lead_time: float            # replenishment lead time (weeks)
    demand_mean_weekly: float   # mean downstream demand per week
    demand_std_weekly: float    # std dev of weekly demand
    holding_cost: float         # $/unit/week at this echelon
    service_level: float        # target fill rate (Clark-Scarf)

def echelon_safety_stock(e: Echelon) -> int:
    z = stats.norm.ppf(e.service_level)
    sigma_lt = e.demand_std_weekly * np.sqrt(e.lead_time)
    return int(np.ceil(z * sigma_lt))

def optimize_network(echelons: List[Echelon]) -> dict:
    """Clark-Scarf decomposition: each echelon optimized independently."""
    results, total_cost = {}, 0
    for e in echelons:
        ss = echelon_safety_stock(e)
        base = int(np.ceil(e.demand_mean_weekly * e.lead_time + ss))
        weekly_cost = ss * e.holding_cost
        results[e.name] = {
            "safety_stock": ss, "base_stock_level": base,
            "weekly_holding_cost": round(weekly_cost, 2),
            "service_level": f"{e.service_level*100:.1f}%"
        }
        total_cost += weekly_cost
    results["total_weekly_ss_cost"] = round(total_cost, 2)
    return results

# Amazon 3-echelon Echo Dot network
network = [
    Echelon("Central Warehouse (ONT8)", 4, 8500, 600, 0.08, 0.99),
    Echelon("Regional DC (BOS1)",        1.5, 2100, 180, 0.12, 0.97),
    Echelon("Fulfillment Center (BOS5)", 2, 1000, 85,  0.16, 0.95),
]
result = optimize_network(network)
for name, data in result.items():
    if name != "total_weekly_ss_cost":
        print(f"{name}: SS={data['safety_stock']}, Base={data['base_stock_level']}")
print(f"Total weekly SS holding cost: ${result['total_weekly_ss_cost']:,.2f}")
# Output:
# Central Warehouse (ONT8): SS=1394, Base=35394
# Regional DC (BOS1): SS=273, Base=3423
# Fulfillment Center (BOS5): SS=197, Base=2197
# Total weekly SS holding cost: $171.98