Single-Echelon Inventory Optimization¶
Stockpyl contains code to solve the following types of single-echelon inventory optimization problems:
The economic order quantity (EOQ) problem (
eoq
module)The newsvendor problem (
newsvendor
module)The (r,Q) optimization problem (
rq
module)The (s,S) optimization problem (
ss
module)The Wagner-Whitin problem (
wagner_whitin
module)Finite-horizon stochastic problems with or without fixed costs (
finite_horizon
module)Plus a number of single-echelon problems with supply uncertainty
Note
The notation and references (equations, sections, examples, etc.) used below refer to Snyder and Shen, Fundamentals of Supply Chain Theory (FoSCT), 2nd edition (2019).
Contents
The EOQ Problem¶
The eoq
module contains code for solving the economic order quantity
(EOQ) problem and some of its variants. The stockpyl.eoq.economic_order_quantity()
function
implements the basic EOQ model; it returns both the optimal order quantity and the corresponding
optimal cost:
>>> from stockpyl.eoq import economic_order_quantity
>>> Q, cost = economic_order_quantity(fixed_cost=8, holding_cost=0.225, demand_rate=1300)
>>> Q
304.0467800264368
>>> cost
68.41052550594829
The module also contains functions for the EOQ with backorders (EOQB) and the economic production quantity (EPQ).
The stockpyl.eoq.joint_replenishment_problem_silver_heuristic()
function implements
Silver’s (1976) heuristic for the joint replenishment problem (JRP):
>>> from stockpyl.eoq import joint_replenishment_problem_silver_heuristic
>>> shared_fixed_cost = 600
>>> individual_fixed_costs = [120, 840, 300]
>>> holding_costs = [160, 20, 50]
>>> demand_rates = [1, 1, 1]
>>> Q, T, m, cost = joint_replenishment_problem_silver_heuristic(shared_fixed_cost, individual_fixed_costs, holding_costs, demand_rates)
>>> Q # order quantities
[3.103164454170876, 9.309493362512628, 3.103164454170876]
>>> T # base cycle time
3.103164454170876
>>> m # order multiples
[1, 3, 1]
>>> cost
837.8544026261366
The Newsvendor Problem¶
The newsvendor
module contains code for solving the newsvendor problem and some of its variants.
The stockpyl.newsvendor.newsvendor_normal()
function
implements the basic newsvendor model for normally distributed demands; it returns both the optimal base-stock level (i.e., order quantity)
and the corresponding expected optimal cost:
>>> from stockpyl.newsvendor import newsvendor_normal
>>> S, cost = newsvendor_normal(holding_cost=0.18, stockout_cost=0.70, demand_mean=50, demand_sd=8)
>>> S
56.60395592743389
>>> cost
1.9976051931766445
If you only want to calculate the expected cost of a given base-stock level, you can either call
stockpyl.newsvendor.newsvendor_normal()
, passing it the optional base_stock_level
parameter, or
call the stockpyl.newsvendor.newsvendor_normal_cost()
function:
>>> from stockpyl.newsvendor import newsvendor_normal, newsvendor_normal_cost
>>> newsvendor_normal(holding_cost=0.18, stockout_cost=0.70, demand_mean=50, demand_sd=8, base_stock_level=53)
(53, 2.223748044859164)
>>> newsvendor_normal_cost(53, holding_cost=0.18, stockout_cost=0.70, demand_mean=50, demand_sd=8)
2.223748044859164
These functions use the version of the newsvendor problem based on holding and stockout (overage and
underage) costs. The stockpyl.newsvendor.newsvendor_normal_explicit()
function uses the “explicit” (or profit-maximization)
form of the newsvendor problem, whose parameters are the selling revenue, purchase cost, and salvage value:
>>> from stockpyl.newsvendor import newsvendor_normal_explicit
>>> newsvendor_normal_explicit(revenue=1.00, purchase_cost=0.30, salvage_value=0.12, demand_mean=50, demand_sd=8)
(56.60395592743389, 33.002394806823354)
The module supports probability distributions other than normal; for example, Poisson:
>>> from stockpyl.newsvendor import newsvendor_poisson
>>> h, p, mu = 0.18, 0.70, 50
>>> newsvendor_poisson(h, p, mu)
(56.0, 1.797235211809178)
You can also solve newsvendor problems for arbitrary continuous distributions specified as scipy.stats.rv_continuous
objects:
>>> from stockpyl.newsvendor import newsvendor_continuous
>>> from scipy.stats import lognorm
>>> from math import exp
>>> demand_distrib = lognorm(0.3, 0, exp(6))
>>> newsvendor_continuous(holding_cost=1, stockout_cost=4, demand_distrib=demand_distrib)
(519.3023987673176, 198.42277610622506)
Or arbitrary discrete distributions specified as either a scipy.stats.rv_discrete
object or
as a dictionary containing the pmf:
>>> from stockpyl.newsvendor import newsvendor_discrete
>>> from scipy.stats import geom
>>> demand_distrib = geom(0.2)
>>> newsvendor_discrete(holding_cost=1, stockout_cost=4, demand_distrib=demand_distrib)
(8.0, 7.194304)
>>> # or ...
>>> d = range(0, 50)
>>> f = [geom.pmf(d_val, 0.2) for d_val in d]
>>> demand_pmf = dict(zip(d, f))
>>> newsvendor_discrete(holding_cost=1, stockout_cost=4, demand_pmf=demand_pmf)
(8, 7.19102133030678)
The \((r,Q)\) Optimization Problem¶
The rq
module contains code for solving the \((r,Q)\) optimization problem, including
a number of approximations.
The stockpyl.rq.r_q_poisson_exact()
function implements Federgruen and Zheng’s (1992)
exact algorithm for the \((r,Q)\) problem with Poisson demands:
>>> from stockpyl.rq import r_q_poisson_exact
>>> r, Q, cost = r_q_poisson_exact(holding_cost=20, stockout_cost=150, fixed_cost=100, demand_mean=1.5, lead_time=2)
>>> r
3
>>> Q
5
>>> cost
107.92358063314975
For normally distributed demands, the module implements the expected-inventory-level (EIL) approximation (Whitin (1953), Hadley and Whitin (1963)), the EOQ with backorders (EOQB) approximation (see Zheng (1992)), the EOQ plus safety stock (EOQ+SS) approximation, and the loss-function approximation (Hadley and Whitin (1963)):
>>> from stockpyl.rq import r_q_eil_approximation, r_q_eoqb_approximation, r_q_eoqss_approximation, r_q_loss_function_approximation
>>> h = 0.225
>>> p = 7.5
>>> K = 8
>>> demand_mean = 1300
>>> demand_sd = 150
>>> L = 1/12
>>> r_eil, Q_eil, _ = r_q_eil_approximation(h, p, K, demand_mean, demand_sd, L)
>>> r_eil, Q_eil
(213.97044213580244, 318.5901810768729)
>>> r_eoqb, Q_eoqb = r_q_eoqb_approximation(h, p, K, demand_mean, demand_sd, L)
>>> r_eoqb, Q_eoqb
(128.63781442427097, 308.5737801203754)
>>> r_eoqss, Q_eoqss = r_q_eoqss_approximation(h, p, K, demand_mean, demand_sd, L)
>>> r_eoqss, Q_eoqss
(190.3369965715624, 304.0467800264368)
>>> r_lf, Q_lf = r_q_loss_function_approximation(h, p, K, demand_mean, demand_sd, L)
>>> r_lf, Q_lf
(126.8670634479628, 328.4491421980451)
We can evaluate these approximate solutions under the exact expected cost function:
>>> from stockpyl.rq import r_q_cost
>>> r_q_cost(r_eil, Q_eil, h, p, K, demand_mean, demand_sd, L)
92.28687665608078
>>> r_q_cost(r_eoqb, Q_eoqb, h, p, K, demand_mean, demand_sd, L)
78.20243187688158
>>> r_q_cost(r_eoqss, Q_eoqss, h, p, K, demand_mean, demand_sd, L)
87.04837003438256
>>> r_q_cost(r_lf, Q_lf, h, p, K, demand_mean, demand_sd, L)
78.07114627035178
The \((s,S)\) Optimization Problem¶
The ss
module contains code for solving the \((s,S)\) optimization problem.
The exact problem with Poisson (or other discrete) demands can be solved using
the stockpyl.ss.s_s_discrete_exact()
function, which implements Zheng and Federgruen’s (1991)
algorithm:
>>> from stockpyl.ss import s_s_discrete_exact
>>> h = 1
>>> p = 4
>>> K = 5
>>> demand_mean = 6
>>> r, Q, cost = s_s_discrete_exact(h, p, K, use_poisson=True, demand_mean=demand_mean)
>>> r
4.0
>>> Q
10.0
>>> cost
8.034111561471642
The stockpyl.ss.s_s_power_approximation()
function implements the power approximation
for normal demands by Ehrhardt and Mosier (1984):
>>> from stockpyl.ss import s_s_power_approximation
>>> h = 0.18
>>> p = 0.70
>>> K = 2.5
>>> demand_mean = 50
>>> demand_sd = 8
>>> r, Q = s_s_power_approximation(h, p, K, demand_mean, demand_sd)
>>> r
40.19461695647407
>>> Q
74.29017010980579
The Wagner-Whitin Problem¶
The wagner_whitin
module contains code for solving the Wagner-Whitin problem using dynamic programming:
>>> from stockpyl.wagner_whitin import wagner_whitin
>>> T = 4
>>> h = 2
>>> K = 500
>>> d = [90, 120, 80, 70]
>>> Q, cost, theta, s = wagner_whitin(T, h, K, d)
>>> Q # Optimal order quantities
[0, 210, 0, 150, 0]
>>> cost # Optimal cost
1380.0
>>> theta # Cost-to-go function
array([ 0., 1380., 940., 640., 500., 0.])
>>> s # Optimal next period to order in
[0, 3, 5, 5, 5]
The cost parameters may also vary from period to period:
>>> h = [5, 1, 1, 2]
>>> K = [300, 500, 300, 500]
>>> Q, cost, _, _ = wagner_whitin(T, h, K, d)
>>> Q
[0, 90, 270, 0, 0]
>>> cost
1020.0
Finite-Horizon Stochastic Problems¶
The finite_horizon
module contains code for solving finite-horizon, stochastic
inventory optimization problems, with or without fixed costs, using dynamic programming (DP).
If the fixed costs are 0, then a base-stock policy is optimal and the results of the demand_pmf indicate the optimal base-stock levels, i.e., order-up-to levels (\(S\)) in each time period:
>>> from stockpyl.finite_horizon import finite_horizon_dp
>>> T = 5
>>> h = 1
>>> p = 20
>>> h_terminal = 1
>>> p_terminal = 20
>>> c = 2
>>> K = 0
>>> mu = 100
>>> sigma = 20
>>> s, S, cost, _, _, _ = finite_horizon_dp(T, h, p, h_terminal, p_terminal, c, K, mu, sigma)
>>> S # Order-up-to levels
[0, 133.0, 133.0, 133.0, 133.0, 126.0]
>>> s # Reorder points equal order-up-to levels in a base-stock policy
[0, 133, 133, 133, 133, 126]
If the fixed costs are non-zero, then an \((s,S)\) policy is optimal, and the results give both the reorder points (\(s\)) and the order-up-to levels (\(S\))
>>> K = 50
>>> s, S, cost, _, _, _ = finite_horizon_dp(T, h, p, h_terminal, p_terminal, c, K, mu, sigma)
>>> s # Reorder points
[0, 110, 110, 110, 110, 111]
>>> S # Order-up-to levels
[0, 133.0, 133.0, 133.0, 133.0, 126.0]