Inventory on a shared shelf

Eight products, each with its own demand distribution and its own cost of running short, share one shelf of 1020 units. For a single product there is a textbook formula; once the products compete for the shelf, there is none. The model below states the expected shortage cost over 500 sampled scenarios the way you would write it in code, max(d − Q, 0) included, and Quicopt solves it as written.

Splitting the shelf in proportion to average demand costs 1347.0 per day. The optimized allocation costs 1014.5, about 24% less, and stays there on scenarios it was not optimized on.

Demand from your sales history

Mean and standard deviation per product, from a daily sales export.

mu_sigma.py
import pandas as pd

sales    = pd.read_csv("sales_daily.csv", index_col=0)   # first column is the date
products = list(sales.columns)
mu       = sales.mean().tolist()            # average daily demand
sigma    = sales.std().tolist()             # how much it varies

The model

inventory_multi.py
import pyomo.environ as pyo
from quicopt import Client
from quicopt.stochastic import Normal, set_distribution, set_scenarios, expectation, maximum

products = ["A", "B", "C", "D", "E", "F", "G", "H"]
mu       = [200.0, 150.0,  90.0, 300.0, 120.0,  60.0, 180.0, 100.0]  # average daily demand
sigma    = [ 20.0,  55.0,  30.0,  30.0,  45.0,  25.0,  18.0,  38.0]  # how much it varies
short    = [  4.0,   4.0,  12.0,   2.0,   9.0,  15.0,   3.0,   7.0]  # cost of one missing unit
hold     = 1.0                                                       # cost of one unit in stock
capacity = 1020.0                                                    # shared shelf space

m = pyo.ConcreteModel()
set_scenarios(m, 500, seed=7)                            # 500 possible days, drawn once

m.Q = pyo.Var(range(8), domain=pyo.NonNegativeIntegers)  # what we order, in whole units
m.d = pyo.Var(range(8))                                  # demand: one random variable per product
for i in range(8):
    set_distribution(m, m.d[i], Normal(mu[i], sigma[i]), name=f"demand_{products[i]}")

m.shelf = pyo.Constraint(expr=sum(m.Q[i] for i in range(8)) <= capacity)   # it has to fit

cost = sum(short[i] * maximum(m.d[i] - m.Q[i], 0.0)
           + hold * maximum(m.Q[i] - m.d[i], 0.0) for i in range(8))
m.cost = pyo.Objective(expr=expectation(cost), sense=pyo.minimize)   # expected cost, all 500 days

result = Client().solve(m)
print(result.display)

What it prints

One run per language. Repeated runs can stop at a slightly different allocation on the same plateau.

▄▀▀▀▄        ▄                     █
█   █ █   █ ▄▄  ▄▀▀▀▄ ▄▀▀▀▄ █▀▀▀▄ ▀█▀▀
█ ▀▄▀ █   █  █  █   ▄ █   █ █▄▄▄▀  █  ▄
 ▀▀ ▀  ▀▀▀▀ ▀▀▀  ▀▀▀   ▀▀▀  █       ▀▀
├── shots
│   ├── 1 · primal   1014.5233   1.0s  ◀ best
│   └── 2 · primal   1014.5233   1.0s
├── status:     heuristic
├── feasible:   true
├── objective:  1014.523297602143
├── x:          x1=193, x2=140, x3=113, x4=73, x5=139, x6=82, …  (8 variables)
└── solve_time: 1.0008 s
planned with the average : 1347.0 EUR/day
optimized                : 1014.5 EUR/day

product   avg demand   by average   optimized
   A          200          170         193
   B          150          128         140
   C           90           76         113
   D          300          255          73
   E          120          102         139
   F           60           51          82
   G          180          153         169
   H          100           85         111
ProductMean demandShortage costProportional splitOptimized
A2004170193
B1504128140
C901276113
D300225573
E1209102139
F60155182
G1803153169
H100785111

Terms you can simply add

Each variant is one more line on the same model, with no reformulation.

A flat fee per stockout event

m.cost = pyo.Objective(
    expr=expectation(cost) + 50.0 * sum(prob(m.d[i] - m.Q[i] >= 0.0) for i in range(8)),
    sense=pyo.minimize)

Cost per day with the fee: 1197.3.

A service guarantee

m.service = pyo.Constraint(
    expr=prob(sum(maximum(m.d[i] - m.Q[i], 0.0) for i in range(8)) <= 350.0) >= 0.90)

Cost per day with the guarantee: 1047.5.