Inventory Optimization & EOQ¶

Ever wondered how Walmart decides when to restock millions of products?
Inventory management is the backbone of supply chain operations. Businesses constantly face two challenges:
1. Overstocking: Too much inventory ties up capital and increases holding costs.
2. Understocking: Too little inventory leads to stockouts, lost sales, and unhappy customers.

To tackle these, companies use Inventory Optimization techniques. In this post, we’ll build a Mixed Integer Programming (MIP) model using Google OR-Tools to:
1. Decide when to order and how much to order
2. Minimize costs while meeting demand forecasts
3. Incorporate Safety Stock, Lead Time, and Economic Order Quantity (EOQ) principles

We’ll use a real dataset from Kaggle: Retail Store Inventory Forecasting Dataset. This dataset contains:
- Date: Daily records
- Store ID, Product ID: Identify SKUs
- Inventory Level: Current stock
- Units Sold: Past sales
- Demand Forecast: Predicted demand
- Price: Item price

We’ll use this to compute:
- Initial inventory
- Safety stock
- Cost parameters

Objective: The objective is to minimize a cost that balances overstocking (holding cost) and understocking (stockout penalty), plus ordering costs.

Understanding Inventory Optimization¶

Inventory optimization is about balancing service level (meeting customer demand) and cost efficiency.

Key concepts¶

Demand Forecast: Predict how much customers will buy.
Safety Stock: Extra inventory to handle uncertainty.
Lead Time: Time between placing an order and receiving it.

Key Costs¶

Holding Cost: Cost of storing inventory.
Stockout Cost: Penalty for not meeting demand.
Ordering Cost: Cost of placing an order.

Why Integer Programming?¶

Decisions like “order or not” are binary.
Quantities are integers (you can’t order 2.5 units).
Constraints like lead time and safety stock make it complex but can be handled by integer variables.

import pandas as pd
import numpy as np
import os
from datetime import timedelta
import matplotlib.pyplot as plt
from ortools.linear_solver import pywraplp
pd.options.display.float_format = '{:,.4f}'.format

Dataset Overview¶

Columns:
1. Date: Daily records.
2. Store ID, Product ID: Identify SKUs.
3. Inventory Level: Current stock.
4. Units Sold: Past sales.
5. Demand Forecast: Predicted demand.
6. Price: Item price.

We can compute initial inventory, forecasted demand, and price-based costs. Historical sales help estimate demand variability for safety stock.

raw_data = pd.read_csv('retail_store_inventory.csv')
# Clean date
raw_data['Date'] = pd.to_datetime(raw_data['Date'])
# Sort by keys
raw_data = raw_data.sort_values(['Store ID','Product ID','Date'])
raw_data

	Date	Store ID	Product ID	Category	Region	Inventory Level	Units Sold	Units Ordered	Demand Forecast	Price	Discount	Weather Condition	Holiday/Promotion	Competitor Pricing	Seasonality
0	2022-01-01	S001	P0001	Groceries	North	231	127	55	135.4700	33.5000	20	Rainy	0	29.6900	Autumn
100	2022-01-02	S001	P0001	Groceries	West	116	81	104	92.9400	27.9500	10	Cloudy	0	30.8900	Spring
200	2022-01-03	S001	P0001	Electronics	West	154	5	189	5.3600	62.7000	20	Rainy	0	58.2200	Winter
300	2022-01-04	S001	P0001	Groceries	South	85	58	193	52.8700	77.8800	15	Cloudy	1	75.9900	Winter
400	2022-01-05	S001	P0001	Groceries	South	238	147	37	150.2700	28.4600	20	Sunny	1	29.4000	Winter
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
72699	2023-12-28	S005	P0020	Groceries	South	198	56	27	50.1800	21.7500	5	Sunny	1	25.2900	Winter
72799	2023-12-29	S005	P0020	Clothing	East	446	268	30	267.5400	85.5800	20	Sunny	1	87.6300	Summer
72899	2023-12-30	S005	P0020	Toys	North	251	149	181	162.9200	79.4800	10	Cloudy	1	82.6900	Autumn
72999	2023-12-31	S005	P0020	Furniture	East	64	40	99	59.6900	90.7900	5	Snowy	1	91.6700	Winter
73099	2024-01-01	S005	P0020	Groceries	East	117	6	165	2.3300	78.3900	20	Rainy	1	79.5200	Spring

73100 rows × 15 columns

In this dataset we are looking at 20 products across 5 stores for two years, from 2022 January 1^st to 2024 January 1^st.

raw_data['Product ID'].unique()

array(['P0001', 'P0002', 'P0003', 'P0004', 'P0005', 'P0006', 'P0007',
       'P0008', 'P0009', 'P0010', 'P0011', 'P0012', 'P0013', 'P0014',
       'P0015', 'P0016', 'P0017', 'P0018', 'P0019', 'P0020'], dtype=object)

raw_data['Store ID'].unique()

array(['S001', 'S002', 'S003', 'S004', 'S005'], dtype=object)

raw_data.Date.min(), raw_data.Date.max()

(Timestamp('2022-01-01 00:00:00'), Timestamp('2024-01-01 00:00:00'))

Calculating Inventory metrics¶

In this section, based on the data we have, we will try to compute the safety stock of the different products and the initial inventory. We will assume that we are on 2023-01-01 or 1^st January 2023, and we will try to optimise inventory for the period after that.

Safety Stock Calculation¶

Safety stock is the stock that should be present in the inventory at all times. This stick halps us if we have sudden increase in demand that the demand forecasting model was unable to explain. It is proportional to the standard deviation of demand (more volatile the demand is, more safety stock is needed) and square root of lead time (more time it takes to recieve a product, more the safety stock). It is given by
$$ SS = z \times \sigma \times \sqrt{L} $$ Where:
- $z$: Service level factor (e.g., 1.65 for 95%).
- $\sigma$: Standard deviation of daily demand.
- $L$: Lead time in days.

# Assumptions

# Assuming we are on January 1st 2023
CUTOFF_DATE = pd.to_datetime('2023-01-01')

# Lead time (days) for all products; can be overridden by mapping later
DEFAULT_LEAD_TIME = 2

# Safety stock service level (z-score): 1.65 ~ 95%, 2.33 ~ 99%
Z_SERVICE = 1.65

# Historical lookback window (days) to estimate demand variability (std dev)
VARIABILITY_LOOKBACK_DAYS = 60

current_levels = raw_data[(raw_data.Date >= CUTOFF_DATE - pd.Timedelta(days=VARIABILITY_LOOKBACK_DAYS-1)) & (raw_data.Date <= CUTOFF_DATE)].\
    groupby(['Product ID', 'Store ID']).\
    aggregate({'Units Sold': ['mean', 'std'],
              'Date': 'max',
              'Units Ordered': 'max',
              'Inventory Level': 'max'}).\
    reset_index()
current_levels.columns = ['product_id', 'store_id', 'mean_sales', 'std_sales', 'date', 'max_units_ordered', 'max_inventory']

# Computing safety stock
current_levels['safety_stock'] = (Z_SERVICE*current_levels['std_sales']*np.sqrt(DEFAULT_LEAD_TIME)).astype(int)

Initial inventory¶

Initial inventory is the inventory as of January 1^st 2023. Based on initial inventory, safety stock and demand, we can predict the future stock.

# Adding current inventory
current_levels = current_levels.merge(raw_data, left_on = ['product_id', 'store_id', 'date'], right_on = ['Product ID', 'Store ID', 'Date'])
current_levels = current_levels[['date', 'product_id', 'store_id', 'safety_stock', 'Inventory Level', 'max_units_ordered', 'max_inventory']]
current_levels

	date	product_id	store_id	safety_stock	Inventory Level	max_units_ordered	max_inventory
0	2023-01-01	P0001	S001	236	425	197	500
1	2023-01-01	P0001	S002	275	104	192	495
2	2023-01-01	P0001	S003	297	219	200	499
3	2023-01-01	P0001	S004	267	117	196	497
4	2023-01-01	P0001	S005	228	430	199	489
...	...	...	...	...	...	...	...
95	2023-01-01	P0020	S001	219	358	200	498
96	2023-01-01	P0020	S002	278	322	199	481
97	2023-01-01	P0020	S003	280	258	199	496
98	2023-01-01	P0020	S004	227	447	199	497
99	2023-01-01	P0020	S005	248	308	197	491

100 rows × 7 columns

EOQ¶

EOQ (Economic Order Quantity) is the optimal order quantity that minimizes total inventory costs. When managing inventory, you face two main costs:
1. Ordering Cost (K): Cost incurred every time you place an order (administration, shipping, setup). (Fewer orders → lower ordering cost)
2. Holding Cost (h): Cost of storing inventory (warehouse space, insurance, depreciation). (Larger orders → higher holding cost)

If we order too frequently, ordering costs skyrocket. If we order too much at once, holding costs explode. EOQ finds the sweet spot where these two costs are balanced.
$$ Q^* = \sqrt{\frac{2DK}{h}} $$
here:
- $Q^*$∗: Optimal order quantity
- D: Annual demand (units/year)
- K: Cost per order (fixed)
- h: Holding cost per unit per year

EOQ assums that:
1. The demand is constant
2. Instant replenishment (lead time=0)
3. No stockouts
4. Single product

Due to these assumptions, it is not possible to use EOQ for complex supply chains directly. Instead we try to optimise the costs and identify when and how much to order.

OR Models¶

OR (Operations Research) models using CP-SAT Solver in OR-Tools are a powerful way to handle combinatorial optimization problems where decisions involve discrete choices, logical conditions, and complex constraints. CP-SAT stands for Constraint Programming with SAT (Boolean Satisfiability). It’s part of Google’s OR-Tools and is designed for:
1. Discrete decision problems (e.g., scheduling, routing, inventory planning).
2. Problems with logical constraints (if-then conditions, Boolean decisions).
3. Large-scale combinatorial optimization.

Unlike linear solvers (which use LP/MIP), CP-SAT uses constraint propagation + SAT solving + branch-and-bound to efficiently explore feasible solutions.

Initialising the model

from ortools.sat.python import cp_model
inventory_model = cp_model.CpModel()

Mathematical Model¶

Decision Variables¶

We have two decisions that we should take, how much to order and when to order.
$order\_rec\_qty_{p,s,t}$: How much order should be recieved (after placing the order) for product p at store q on day t.
$order\_bool_{p,s,t}$: Whether to place an order (binary).

future_data = raw_data[(raw_data['Store ID'] == 'S001') & (raw_data['Product ID'] == 'P0001') & (raw_data.Date>=CUTOFF_DATE)]
future_data = future_data[['Date', 'Store ID', 'Product ID', 'Demand Forecast']]
future_data.columns = ['date', 'store_id', 'product_id', 'demand']
future_data = future_data.merge(current_levels, on=['date', 'store_id', 'product_id'], how='left')

# Assuming constant safety stock
future_data['safety_stock'] = future_data.groupby(['store_id', 'product_id'])['safety_stock'].ffill().astype(int)

# Assuming max units to order doesnt change
future_data['max_units_ordered'] = future_data.groupby(['store_id', 'product_id'])['max_units_ordered'].ffill().astype(int)

# Assuming max inventory is the capasity of the product inventory in the store
future_data['max_inventory'] = future_data.groupby(['store_id', 'product_id'])['max_inventory'].ffill().astype(int)

# Converting demand to integer
future_data['demand'] = future_data.demand.astype(int)

# Creating one boolean decision variable for whether to order product-store at time t to order
future_data['order_bool'] = None
for i in range(len(future_data)):
    date, store, prod = (future_data.date.dt.date[i], future_data.store_id[i], future_data.product_id[i])
    future_data.loc[i, 'order_bool'] = inventory_model.NewBoolVar('order_bool_%s,%s,%s' % (date, store, prod))
    if i<5:
        print('Creating the boolean variable ', future_data.loc[i, 'order_bool'], 
          'representing the if we should order ', date, store, prod)

Creating the boolean variable  order_bool_2023-01-01,S001,P0001 representing the if we should order  2023-01-01 S001 P0001
Creating the boolean variable  order_bool_2023-01-02,S001,P0001 representing the if we should order  2023-01-02 S001 P0001
Creating the boolean variable  order_bool_2023-01-03,S001,P0001 representing the if we should order  2023-01-03 S001 P0001
Creating the boolean variable  order_bool_2023-01-04,S001,P0001 representing the if we should order  2023-01-04 S001 P0001
Creating the boolean variable  order_bool_2023-01-05,S001,P0001 representing the if we should order  2023-01-05 S001 P0001

# Creating one integer decision variable for how much to order for each product-store at time t to order
future_data['order_rec_qty'] = None
for i in range(len(future_data)):
    date, store, prod = (future_data.date.dt.date[i], future_data.store_id[i], future_data.product_id[i])
    future_data.loc[i, 'order_rec_qty'] = inventory_model.NewIntVar(0, future_data.max_units_ordered[i]*10, 'order_rec_qty_%s,%s,%s' % (date, store, prod))
    if i<5:
        print('Creating the integer variable ', future_data.loc[i, 'order_rec_qty'], 
          'representing how much we should order ', date, store, prod)

Creating the integer variable  order_rec_qty_2023-01-01,S001,P0001 representing how much we should order  2023-01-01 S001 P0001
Creating the integer variable  order_rec_qty_2023-01-02,S001,P0001 representing how much we should order  2023-01-02 S001 P0001
Creating the integer variable  order_rec_qty_2023-01-03,S001,P0001 representing how much we should order  2023-01-03 S001 P0001
Creating the integer variable  order_rec_qty_2023-01-04,S001,P0001 representing how much we should order  2023-01-04 S001 P0001
Creating the integer variable  order_rec_qty_2023-01-05,S001,P0001 representing how much we should order  2023-01-05 S001 P0001

$inventory_{p,s,t}$: Inventory at end of day

# Creating one integer decision variable for inventory
future_data['inventory'] = None
for i in range(len(future_data)):
    date, store, prod = (future_data.date.dt.date[i], future_data.store_id[i], future_data.product_id[i])
    future_data.loc[i, 'inventory'] = inventory_model.NewIntVar(0, future_data.max_inventory[i]*2, 'inventory_%s,%s,%s' % (date, store, prod))
    if i<5:
        print('Creating the integer variable ', future_data.loc[i, 'inventory'], 
          'representing how much inventory is present ', date, store, prod)

Creating the integer variable  inventory_2023-01-01,S001,P0001 representing how much inventory is present  2023-01-01 S001 P0001
Creating the integer variable  inventory_2023-01-02,S001,P0001 representing how much inventory is present  2023-01-02 S001 P0001
Creating the integer variable  inventory_2023-01-03,S001,P0001 representing how much inventory is present  2023-01-03 S001 P0001
Creating the integer variable  inventory_2023-01-04,S001,P0001 representing how much inventory is present  2023-01-04 S001 P0001
Creating the integer variable  inventory_2023-01-05,S001,P0001 representing how much inventory is present  2023-01-05 S001 P0001

$stockout_{p,s,t}$ stockout during the day

future_data['stockout'] = None
for i in range(len(future_data)):
    date, store, prod = (future_data.date.dt.date[i], future_data.store_id[i], future_data.product_id[i])
    future_data.loc[i, 'stockout'] = inventory_model.NewIntVar(0, future_data.max_inventory[i]*2, 'stockout_%s,%s,%s' % (date, store, prod))

future_data

	date	store_id	product_id	demand	safety_stock	Inventory Level	max_units_ordered	max_inventory	order_bool	order_rec_qty	inventory	stockout
0	2023-01-01	S001	P0001	194	236	425.0000	197	500	order_bool_2023-01-01,S001,P0001	order_rec_qty_2023-01-01,S001,P0001	inventory_2023-01-01,S001,P0001	stockout_2023-01-01,S001,P0001
1	2023-01-02	S001	P0001	98	236	NaN	197	500	order_bool_2023-01-02,S001,P0001	order_rec_qty_2023-01-02,S001,P0001	inventory_2023-01-02,S001,P0001	stockout_2023-01-02,S001,P0001
2	2023-01-03	S001	P0001	248	236	NaN	197	500	order_bool_2023-01-03,S001,P0001	order_rec_qty_2023-01-03,S001,P0001	inventory_2023-01-03,S001,P0001	stockout_2023-01-03,S001,P0001
3	2023-01-04	S001	P0001	261	236	NaN	197	500	order_bool_2023-01-04,S001,P0001	order_rec_qty_2023-01-04,S001,P0001	inventory_2023-01-04,S001,P0001	stockout_2023-01-04,S001,P0001
4	2023-01-05	S001	P0001	152	236	NaN	197	500	order_bool_2023-01-05,S001,P0001	order_rec_qty_2023-01-05,S001,P0001	inventory_2023-01-05,S001,P0001	stockout_2023-01-05,S001,P0001
...	...	...	...	...	...	...	...	...	...	...	...	...
361	2023-12-28	S001	P0001	78	236	NaN	197	500	order_bool_2023-12-28,S001,P0001	order_rec_qty_2023-12-28,S001,P0001	inventory_2023-12-28,S001,P0001	stockout_2023-12-28,S001,P0001
362	2023-12-29	S001	P0001	184	236	NaN	197	500	order_bool_2023-12-29,S001,P0001	order_rec_qty_2023-12-29,S001,P0001	inventory_2023-12-29,S001,P0001	stockout_2023-12-29,S001,P0001
363	2023-12-30	S001	P0001	20	236	NaN	197	500	order_bool_2023-12-30,S001,P0001	order_rec_qty_2023-12-30,S001,P0001	inventory_2023-12-30,S001,P0001	stockout_2023-12-30,S001,P0001
364	2023-12-31	S001	P0001	42	236	NaN	197	500	order_bool_2023-12-31,S001,P0001	order_rec_qty_2023-12-31,S001,P0001	inventory_2023-12-31,S001,P0001	stockout_2023-12-31,S001,P0001
365	2024-01-01	S001	P0001	53	236	NaN	197	500	order_bool_2024-01-01,S001,P0001	order_rec_qty_2024-01-01,S001,P0001	inventory_2024-01-01,S001,P0001	stockout_2024-01-01,S001,P0001

366 rows × 12 columns

Constraints¶

Inventory constraint¶

inventory at the end of the day depends on
- Previous day inventory
- - current day demand
- + order recieved on the day
- + stockout value of the day (variable included to handle stockouts)

\[ I_{p,s,t} = I_{p,s,t-1} - demand_{p,s,t} + order\_value\_var_{p,s,t} + stockout_{p,s,t}\]

# At every store and every product, inventory equation should be held
for store in future_data.store_id.unique():
    for product in future_data.product_id.unique():
        store_prod_data = future_data[(future_data.store_id == store) & (future_data.product_id == product)]
        store_prod_data = store_prod_data.sort_values('date').reset_index()
        for time in range(len(store_prod_data)):#.date.unique(): # 365 days
            if(time==0):
                lhs = store_prod_data.loc[time, 'inventory']
                rhs = int(store_prod_data.loc[time, 'Inventory Level'])
                inventory_model.Add(lhs == rhs)
            else:
                lhs = store_prod_data.loc[time, 'inventory']
                rhs = store_prod_data.loc[time-1, 'inventory'] - store_prod_data.loc[time, 'demand'] + \
                      store_prod_data.loc[time, 'order_rec_qty'] + store_prod_data.loc[time, 'stockout']
                inventory_model.Add(lhs == rhs)

Order constraint¶

The boolean order_bool should be true when $order\_rec\_qty>0$. This can be represented as:

\[ order\_rec\_qty_{p,s,t} \leq M\times order\_bool_{p,s,t} \]

\[ order\_rec\_qty_{p,s,t} \geq 0.1-M\times (1-order\_bool_{p,s,t}) \]

M = 100000

for i in range(len(future_data)):
    inventory_model.Add(future_data.loc[i, 'order_rec_qty'] <= M*future_data.loc[i, 'order_bool'])
    inventory_model.Add(future_data.loc[i, 'order_rec_qty'] >= 1 - M*(1-future_data.loc[i, 'order_bool']))

Stockout constraints¶

Stockout happens only when inventory is zero at time t. For the constraint $stockout_{p,s,t}>0$ when $inventory_{p,s,t}=0$ and $stockout_{p,s,t}=0$ when $inventory_{p,s,t}>0$ we have to create a boolean $stock\_bool_{p,s,t}$ such that

\[ stockout_{p,s,t} \leq M\times stock\_bool_{p,s,t} \]

\[ inventory_{p,s,t} \leq M\times (1-stock\_bool_{p,s,t}) \]

Creating the boolean variable $stock\_bool_{p,s,t}$ which represents if there is a stock out.

future_data['stock_bool'] = None
for i in range(len(future_data)):
    future_data.loc[i, 'stock_bool'] = inventory_model.NewBoolVar('stock_bool_%s,%s,%s' % (date, store, prod))
    inventory_model.Add(future_data.loc[i, 'stockout'] <= M*future_data.loc[i, 'stock_bool'])
    inventory_model.Add(future_data.loc[i, 'inventory'] <= M*(1-future_data.loc[i, 'stock_bool']))

Safety stock constraints¶

We need to create a soft constraint that inventory should be above the safety stock as much as possible. For this we are creating a slack of how much inventory is above safety stock. We can give less value to the weight of this boolean while minimizing to make this a soft constraint.

\[ safety\_stock\_slack_{p,s,t} >= safety\_stock_{p,s,t} - inventory_{p,s,t} \]

future_data['safety_stock_slack'] = None
for i in range(len(future_data)):
    future_data.loc[i, 'safety_stock_slack'] = inventory_model.NewIntVar(0, future_data.max_inventory[i]*2, 'safety_stock_slack_%s,%s,%s' % (date, store, prod))
    inventory_model.Add(future_data.loc[i, 'safety_stock_slack'] >= future_data.loc[i, 'safety_stock'] - future_data.loc[i, 'inventory'])

Maximum inventory constraint¶

Maximum inventory should be less than the capacity.
$$ \sum_{s}{inventory_{p,s,t}} >= max_inventory_{s} $$

for store in future_data.store_id.unique():
    store_data = future_data[(future_data.store_id == store)]

    rhs = raw_data[raw_data['Store ID'] == store].groupby('Date')['Inventory Level'].sum().max()
    for t in future_data.date.unique():
        lhs = store_data.loc[store_data.date == t, 'inventory'].sum()
        inventory_model.Add(lhs <= rhs)

Objective¶

The objective is to minimise costs. Costs are of four types:
- Holding cost
- stockout penalty
- ordering cost
- penalty for going below safety stock

HOLDING_COST = 0.1
STOCKOUT_PENALTY = 10000
ORDERING_COST = 10000
SAFETY_STOCK_PENALTY = 100

total_cost = sum(future_data['inventory']*(HOLDING_COST)) + sum(future_data['stock_bool']*STOCKOUT_PENALTY) + \
    sum(future_data['order_bool']*ORDERING_COST) + sum(future_data['safety_stock_slack']*SAFETY_STOCK_PENALTY)
#                  /future_data.max_inventory)

# The objective is to minimise cost
inventory_model.Minimize(total_cost)

# Solving the problem
solver = cp_model.CpSolver()
solution_printer = cp_model.ObjectiveSolutionPrinter()
solver.parameters.max_time_in_seconds = 60.0 # ends the solver after 60 seconds
status = solver.SolveWithSolutionCallback(inventory_model, solution_printer)

Solution 0, time = 0.22 s, objective = 12087082
Solution 1, time = 0.25 s, objective = 12008958
Solution 2, time = 0.27 s, objective = 11930183
Solution 3, time = 0.27 s, objective = 2974744
Solution 4, time = 0.34 s, objective = 2456465
Solution 5, time = 0.61 s, objective = 2456446
Solution 6, time = 0.84 s, objective = 2456347
Solution 7, time = 0.92 s, objective = 763137
Solution 8, time = 1.25 s, objective = 763072
Solution 9, time = 2.15 s, objective = 747383
Solution 10, time = 2.31 s, objective = 686536
Solution 11, time = 3.67 s, objective = 655250
Solution 12, time = 4.79 s, objective = 641887
Solution 13, time = 5.35 s, objective = 639389
Solution 14, time = 6.25 s, objective = 639376
Solution 15, time = 7.16 s, objective = 638026
Solution 16, time = 10.30 s, objective = 633782
Solution 17, time = 12.40 s, objective = 617881
Solution 18, time = 14.50 s, objective = 617503
Solution 19, time = 16.86 s, objective = 615899
Solution 20, time = 17.47 s, objective = 615501
Solution 21, time = 20.66 s, objective = 612485
Solution 22, time = 24.29 s, objective = 604623

cp_model.INFEASIBLE == status

False

cp_model.OPTIMAL == status

False

The solution is optimal. The solution is

future_data_solver = future_data.copy()
for i in range(len(future_data_solver)):
    date, store, prod = (future_data.date.dt.date[i], future_data.store_id[i], future_data.product_id[i])
    future_data_solver.loc[i, 'order_bool'] = solver.Value(future_data.loc[i, 'order_bool'])
    future_data_solver.loc[i, 'order_rec_qty'] = solver.Value(future_data.loc[i, 'order_rec_qty'])
    future_data_solver.loc[i, 'inventory'] = solver.Value(future_data.loc[i, 'inventory'])
    future_data_solver.loc[i, 'stockout'] = solver.Value(future_data.loc[i, 'stockout'])
    future_data_solver.loc[i, 'safety_stock_slack'] = solver.Value(future_data.loc[i, 'safety_stock_slack'])
    future_data_solver.loc[i, 'stock_bool'] = solver.Value(future_data.loc[i, 'stock_bool'])

Visualization¶

After solving, plot:
- Inventory vs Demand vs Orders
- Highlight safety stock line
The inventory for the product P001 and store S001 is

p_11 = future_data_solver[(future_data_solver.product_id == 'P0001') & (future_data_solver.store_id == 'S001')]
fs, axs = plt.subplots(3, figsize=(20,30), gridspec_kw={'height_ratios': [2, 1, 1]}, constrained_layout=True)
p_11.plot(x='date', y='inventory', ax=axs[0])
p_11.plot(x='date', y='order_rec_qty', ax=axs[1])
p_11.plot(x='date', y='stockout', ax=axs[2])

axs[0].set_ylabel('Inventory', fontsize=15); axs[0].set_xlabel('', fontsize=25); 
axs[0].set_title('Inventory of P0001 in S001', fontsize=15)
axs[1].set_ylabel('Order recieved Quantity', fontsize=15); axs[2].set_xlabel('', fontsize=15);
axs[2].set_ylabel('Stockouts', fontsize=15); axs[2].set_xlabel('Time', fontsize=15);

plt.show();

png

The chart clearly shows inventory levels fluctuating over time. Notice how the system strives to keep inventory above the safety stock threshold to maintain service levels. When demand consumes stock, inventory gradually declines until a replenishment order is triggered, restoring levels back above the safety stock line. Finally the answer to when to order and how much to order is

future_data_solver['when_to_order'] = future_data_solver.groupby(['store_id', 'product_id']).order_bool.shift(-DEFAULT_LEAD_TIME)
future_data_solver['how_much_to_order'] = future_data_solver.groupby(['store_id', 'product_id']).order_rec_qty.shift(-DEFAULT_LEAD_TIME)
future_data_solver.loc[future_data_solver.when_to_order==1, ['date', 'how_much_to_order']]

	date	how_much_to_order
0	2023-01-01	789
8	2023-01-09	833
13	2023-01-14	883
21	2023-01-22	1084
28	2023-01-29	704
34	2023-02-04	906
39	2023-02-09	1122
43	2023-02-13	1070
51	2023-02-21	999
58	2023-02-28	777
62	2023-03-04	1006
68	2023-03-10	990
78	2023-03-20	916
82	2023-03-24	962
87	2023-03-29	806
90	2023-04-01	939
99	2023-04-10	996
106	2023-04-17	625
113	2023-04-24	770
119	2023-04-30	706
127	2023-05-08	718
133	2023-05-14	809
137	2023-05-18	909
146	2023-05-27	815
154	2023-06-04	930
162	2023-06-12	999
167	2023-06-17	864
174	2023-06-24	1044
180	2023-06-30	924
190	2023-07-10	828
196	2023-07-16	907
200	2023-07-20	908
206	2023-07-26	1171
212	2023-08-01	1057
218	2023-08-07	1016
224	2023-08-13	1049
231	2023-08-20	935
237	2023-08-26	936
247	2023-09-05	1141
258	2023-09-16	1029
263	2023-09-21	1053
269	2023-09-27	822
276	2023-10-04	1161
287	2023-10-15	878
292	2023-10-20	828
298	2023-10-26	843
301	2023-10-29	980
309	2023-11-06	1070
318	2023-11-15	760
321	2023-11-18	799
326	2023-11-23	859
331	2023-11-28	994
339	2023-12-06	897
343	2023-12-10	842
350	2023-12-17	856
355	2023-12-22	1081

References¶

Kaggle dataset: Retail Store Inventory Forecasting Dataset (link).
OR-Tools MIP: https://developers.google.com/optimization.