学習目標
新聞売り子問題のコスト関数を定式化できる.
コストの期待値を最小化することで,最適発注量を導出できる.
臨界率とサービスレベルの関係を理解する.
需要が正規分布に従う場合の最適発注量を求めることができる.
*販売価格や残存価額を含めた再定式化を理解する.
新聞売り子問題とは
これまで紹介した EOQ 在庫モデルは,需要が決定論的であると仮定していた.ここからは,需要が確率的であると仮定した在庫モデルを紹介する.
新聞売り子問題(Newsvendor Problem)は,古典的な確率的在庫モデルの一つである.新聞は次の日には売れなくなるため,新聞売り子問題は最も単純な perishable 在庫モデルとして知られている.新聞に限らず,食品や花などの生鮮品の在庫管理にも応用される.
新聞売り子が新聞を仕入れ,販売する問題を考える.新聞 1 部の欠品費用を p ,保管費用を h とする.新聞売子問題において,p を在庫不足費用 (underage cost),h を在庫超過費用 (overage cost)とも呼ぶ.p \gt 0 ,h \gt 0 とし,初期在庫は 0 とする.
新聞の需要 D を確率変数とするとき,新聞売り子はどれだけの新聞を発注すればよいかという問題である.
練習 17.1 (思考実験) 需要 D が 一様分布 U(100,300) に従うとする.以下のそれぞれの場合,新聞売り子がどれだけの新聞を発注するか考えよ.
在庫超過費用 h = 1000 ,在庫不足費用 p = 0.1
在庫超過費用 h = 0.1 ,在庫不足費用 p = 1000
在庫超過費用 h = 10 ,在庫不足費用 p = 10
記号
h
1 個当たりの在庫超過費用
p
1 個当たりの在庫不足費用
D
需要(確率変数)
d
需要の観測値
f_D(d)
需要 D の確率密度関数
F_D(d)
需要 D の累積分布関数
S
発注量
g(S, d)
発注量 S ,需要の観測値 d に対するコスト
g(S)
g(S, d) の期待値,g(S) = \mathbb{E}[g(S, D)]
定式化
新聞売り子が S 部の新聞を仕入れ,需要 D が d であったとする.このとき,新聞売り子のコスト g(S, d) は以下のように表される.
g(S, d) = h (S - d)^+ + p (d - S)^+
\tag{17.1}
例 17.1 発注量 S = 100 ,需要 D の観測値 d = 120 ,在庫超過費用が h = 10 ,在庫不足費用が p = 5 のとき,コスト g(S, d) は以下のように求められる.
g(100, 120) = 10 (100 - 120)^+ + 5 (120 - 100)^+ = 100
観測値 d が分かれば,コスト g(S, d) を計算できる.しかし,需要 D が確率変数であることを思い出そう.ここからは,需要 D が連続型確率変数であると仮定する.
需要 D の確率密度関数を f_D(d) ,累積分布関数を F_D(d) とする.発注量を S としたとき,需要 D \leq S である確率は,累積分布関数 F_D(S) で与えられる.
\mathbb{P}(D \leq S) = F_D(S) = \int_{0}^{S} f_D(d) dd
\tag{17.2}
例 17.2 需要 D が一様分布 U(100,300) に従うとき,確率密度関数 f_D(250) と累積分布関数 F_D(250) は以下のように表される.
コード
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import uniform
# 一様分布のパラメータ
a = 100
b = 300
# 需要の範囲
d = np.linspace(a - 50 , b + 50 , 1000 )
# plot pdf with shaded area
plt.figure(figsize= (4 , 3.2 ))
plt.plot(d, uniform.pdf(d, a, b - a), color= "#144f78" )
plt.fill_between(
d,
0 ,
uniform.pdf(d, a, b - a),
where= (d <= 250 ),
color= "#144f78" ,
alpha= 0.18 ,
label= "F(250)" ,
)
plt.xlabel(r" $ D $ " )
plt.ylabel(r" $ f_D ( d ) $ " )
plt.xlim(a - 50 , b + 50 )
plt.ylim(0 , 0.006 )
plt.show()
# plot cdf
plt.figure(figsize= (4 , 3.2 ))
plt.plot(d, uniform.cdf(d, a, b - a), color= "#144f78" )
plt.axvline(250 , color= "#b03535" , linestyle= "--" , label= "d=250" )
plt.xlabel(r" $ D $ " )
plt.ylabel(r" $ F_D ( d ) $ " )
plt.xlim(a - 50 , b + 50 )
plt.ylim(0 , 1.1 )
plt.show()
需要 D \leq 250 である確率は,累積分布関数 F_D(250) で与えられる.
\mathbb{P}(D \leq 250) = F_D(250) = \int_{100}^{250} f_D(d) dd = 0.75
需要 D > 250 である確率は 1 - F_D(250) で与えられる.
\mathbb{P}(D > 250) = 1 - F_D(250) = 1 - \int_{100}^{250} f_D(d) dd = 0.25
D が確率変数であるため,コストも確率変数となる.新聞売り子問題の目的は,コストの期待値を最小化することである.
ここで,コストの期待値を g(S) = \mathbb{E}[g(S, D)] とする.g(S, D) は 式 17.1 で与えられるため,コストの期待値 g(S) は以下のように表される.
\begin{align}
g(S) &= \mathbb{E}[g(S, D)] \\
&= \mathbb{E}[h (S - D)^+ + p (D - S)^+] \\
&= \mathbb{E}[h (S - D)^+] + \mathbb{E}[p (D - S)^+] \\
&= h \mathbb{E}[(S - D)^+] + p \mathbb{E}[(D - S)^+] \\
\end{align}
この式により,コストの期待値 g(S) は,(S - D)^+ の期待値かける在庫超過費用 h と,(D - S)^+ の期待値かける在庫不足費用 p の和であることが分かる.
X と Y を確率変数,a と b を定数とする.このとき,
\mathbb{E}[X + Y] = \mathbb{E}[X] + \mathbb{E}[Y]
\mathbb{E}[aX + b] = a\mathbb{E}[X] + b
が成り立つ.
需要 D が連続型確率変数であるため,\mathbb{E}[(S - D)^+] と \mathbb{E}[(D - S)^+] は以下のように表される. \begin{align}
\mathbb{E}[(S - D)^+]
&= \int_{0}^{\infty} (S - d)^+ f_D(d) dd = \int_{0}^{S} (S - d) f_D(d) dd \\
\mathbb{E}[(D - S)^+]
&= \int_{0}^{\infty} (d - S)^+ f_D(d) dd = \int_{S}^{\infty} (d - S) f_D(d) dd
\end{align}
X を連続型確率変数,f_X(x) を X の確率密度関数とする.このとき,X の期待値 \mathbb{E}[X] は以下のように表される.
\mathbb{E}[X] = \int_{-\infty}^{\infty} x f_X(x) dx
\mathbb{E}[g(X)] は以下のように表される.
\mathbb{E}[g(X)] = \int_{-\infty}^{\infty} g(x) f_X(x) dx
したがって,コストの期待値 g(S) は以下のように表される.
\begin{align}
g(S) &= \mathbb{E}[g(S, D)] \\
&= h \mathbb{E}[(S - D)^+] + p \mathbb{E}[(D - S)^+] \\
&= h \int_{0}^{S} (S - d) f_D(d) dd + p \int_{S}^{\infty} (d - S) f_D(d) dd
\end{align}
例 17.3 在庫超過費用が h = 10 ,在庫不足費用が p = 5 のとき,発注量 S = 150 に対するコストの期待値 g(150) は以下のように求められる.
g(150) = 10 \int_{0}^{150} (150 - d) f_D(d) dd + 5 \int_{150}^{\infty} (d - 150) f_D(d) dd
最適化
以上の議論から,新聞売り子問題は,コストの期待値 g(S) を最小化する発注量 S を求める問題に帰着される.手順は以下の通りである.
g(S) の1階微分を求める.
dg(S)/dS = 0 を解く.
2階微分を求め,g(S) が凸関数であることを確認する.
式 17.2 を用い,g(S) の1階微分は以下のように求める.
\begin{align}
\frac{dg(S)}{dS} &= h \int_{0}^{S} f_D(d) dd - p \int_{S}^{\infty} f_D(d) dd \\
&= h F_D(S) - p (1 - F_D(S)) \\
\end{align}
よって,dg(S)/dS = 0 から,
\begin{align}
h F_D(S) - p (1 - F_D(S)) &= 0 \\
F_D(S) &= \frac{p}{h + p}
\end{align}
になる.2 階微分は
\begin{align}
\frac{d^2g(S)}{dS^2} &= h f_D(S) + p f_D(S) \\
&= (h + p) f_D(S) > 0
\end{align}
である.したがって,コスト関数 g(S) は凸関数であり,1 階微分が 0 になる点は最小値を与える.
p > 0 ,h > 0 ,f_D(S) > 0 より,(h + p) f_D(S) > 0 が成り立つ.
コスト関数 g(S) を最小化するための最適発注量 S^* は
S^* = F_D^{-1}\left(\frac{p}{h + p}\right)
となる.ここで,F_D^{-1} は需要 D の累積分布関数の逆関数である.
ある関数 y = f(x) に対し,次の条件を満たす関数 x = f^{-1}(y) を f(x) の逆関数と呼ぶ.
f(f^{-1}(y)) = y, \quad f^{-1}(f(x)) = x
定理 17.1 新聞売り子問題における最適発注量 S^* は,
S^* = F_D^{-1}\left(\frac{p}{h + p}\right)
で与えられる.
正規分布の場合
需要 D が正規分布 N(\mu, \sigma^2) に従うとき,新聞売り子問題の最適発注量 S^* は以下の式を満たす.
F_D(S^*) = \Phi\left(\frac{S^* - \mu}{\sigma}\right) = \frac{p}{h + p}
ここで,
z = \Phi^{-1}\left(\frac{p}{h + p}\right)
とおくと,
S^* = \sigma z + \mu
で与えられる.標準正規分布表を用いて z を調べることができる.
平均 0,分散 1 の正規分布を標準正規分布 (standard normal distribution)と呼ぶ.標準正規分布に従う確率変数を Y とすると,Y \sim N(0, 1) と表される.標準正規分布の確率密度関数を \phi(z) ,累積分布関数を \Phi(z) と表す.
与えられた X \sim N(\mu, \sigma^2) の累積分布関数 F_X(x) の値を求めるには,以下のように変換する.
\begin{align}
F_X(x) &= \mathbb{P}(X \leq x) \\
&= \mathbb{P} \left( \frac{X - \mu}{\sigma} \leq \frac{x - \mu}{\sigma} \right) \\\\
&= \mathbb{P}\left( Y \leq \frac{x - \mu}{\sigma} \right) \\\\
&= \Phi \left( \frac{x - \mu}{\sigma} \right)
\end{align}
与えられた z に対し,\Phi(z) の値は標準正規分布表 (standard normal table)を用いて調べることができる.
臨界率
以上の議論から,最適発注量 S^* は以下の式を満たす.
F_D(S) = \mathbb{P}(D \leq S) = \frac{p}{h + p}
ここで,F_D(S) = \mathbb{P}(D \leq S) は需要 D が発注量 S 以下である確率を表す.言い換えると,欠品が発生しない確率を表す.この確率のことはサービスレベル (service level)と呼ぶ.定理 17.1 は,サービスレベルを p/(h+p) に等しくする発注量 S が最適であることを示している.
この p/(h+p) は臨界率 (critical ratio)と呼ばれる.
例題
例 17.4 需要 D が正規分布 N(100, 25) に従う新聞売り子問題を考える.在庫超過費用が h = 10 ,在庫不足費用が p = 40 のとき,最適発注量 S^* を求める.
定理 17.1 より,S^* は以下の式で与えられる.
S^* = F_D^{-1}\left(\frac{p}{h + p}\right) = F_D^{-1}\left(\frac{40}{10 + 40}\right) = F_D^{-1}\left(0.8\right)
言い換えると,F_D(S^*) = 0.8 を満たす S^* を求めればよい.これを図で表すと,面積が 0.8 になるような S^* を求めることに相当する.
コード
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
# 正規分布のパラメータ
mu = 100
sigma = 5
# 需要の範囲
d = np.linspace(mu - 4 * sigma, mu + 4 * sigma, 1000 )
S_star = norm.ppf(0.8 , loc= mu, scale= sigma)
# plot the normal distribution
plt.plot(d, norm.pdf(d, mu, sigma), color= "#144f78" )
plt.fill_between(
d,
0 ,
norm.pdf(d, mu, sigma),
where= (d <= S_star),
color= "#144f78" ,
alpha= 0.18 ,
label= "F(S*)" ,
)
plt.xlabel("Demand" )
plt.ylabel("Probability Density" )
plt.axvline(mu, color= "#5b6b7c" , linestyle= "--" , label= "mu" )
plt.axvline(S_star, color= "#b03535" , linestyle= "--" , label= "S*" )
plt.ylim(0 , 0.09 )
plt.legend()
plt.show()
下では,Excel,Python,標準正規分布表を用いて S^* を求める方法を示す.
標準正規分布表
標準正規分布表 を調べると,\Phi(z) = 0.8 のとき,z は約 0.84 である.したがって,S^* は
S^* = z \sigma + \mu \approx 0.84 \times 5 + 100 = 104.2
となる.
Excel
Excel では,NORM.INV(確率,平均,標準偏差) 関数を用いて F_D^{-1}(0.8) を求めることができる.
=NORM.INV(0.8, 100, 5)
Python
Python では,SciPy ライブラリの ppf() 関数を用いて,逆関数を求めることができる.
from scipy.stats import norm
# 正規分布のパラメータ
mu = 100
sigma = 5
# 在庫超過費用と在庫不足費用
h = 10
p = 40
# 臨界率
critical_ratio = p / (h + p)
# 最適発注量
S_star = norm.ppf(critical_ratio, loc= mu, scale= sigma)
print (f"Optimal order quantity S*: { S_star:.2f} " )
Optimal order quantity S*: 104.21
再定式化
新聞 1 部の仕入れ価格を c ,販売価格を r ,残存価額を v ,欠品費用を p ,保管費用を h とする.ここで,以下の条件を満たすとする.
r > c . 販売価格は仕入れ価格より高い.
r > v . 販売価格は残存価値より高い.
新聞売り子が S 部の新聞を仕入れ,需要 D が d であったとする.このとき,新聞売り子の利益 \pi(S, d) は以下のように表される.
コスト関数
\begin{align}
\pi(S, d) & = r \min\{d, S\} - cS + v \max\{0, S - d\} \\
& \quad - h \max\{0, S - d\} - p \max\{0, d - S\}
\end{align}
第一項は販売利益,第二項は仕入れコスト,第三項は残存価値,第四項は保管コスト,第五項は欠品コストである.
以下は,\max\{0, x\} = x^+ を用いて書き換えた形である.整理すると,利益は以下のように表される.
\pi(S, d) = r \min\{d, S\} - cS + (v - h) (S - d)^+ - p (d - S)^+
第一項を以下のように書き換えることができる.
r \min\{d, S\} = r d - r (d - S)^+
d < S の場合,r \min\{d, S\} = r d となる.
d >= S の場合,r \min\{d, S\} = r d - r(d - S) = r S となる.
したがって,利益は以下のように書き換えられる.
\pi(S, d) = r d - cS + (v- h) (S - d)^+ - (p + r) (d - S)^+
利益の最大化は,コストの最小化に帰着される.したがって,コスト関数 g(S, d) = -\pi(S, d) は以下のように表される.
g(S, d) = cS - rd + (h - v) (S - d)^+ + (p + r) (d - S)^+
D が確率変数であるため,コストの期待値 g(S) は以下のように表される.
\begin{align}
g(S) &= \mathbb{E}[g(S, D)] \\
&= \int_{0}^{\infty} g(S, d) f_D(d) dd \\
&= cS - r\mathbb{E}[D] + (h - v) \mathbb{E}[(S - D)^+] + (p + r) \mathbb{E}[(D - S)^+] \\
&= c S - r \mu + (h - v) \int_{0}^{\infty}(S - d)^+ f_D(d) dd + (p + r) \int_{0}^{\infty} (d - S)^+ f_D(d) dd \\
&= c S - r \mu + (h - v) \int_{0}^{S} (S - d) f_D(d) dd + (p + r) \int_{S}^{\infty} (d - S) f_D(d) dd
\end{align}
最適化
g(S) の1階微分は以下のように求める.
\begin{align}
\frac{dg(S)}{dS} &= c + (h - v) F_D(S) - (p + r) (1 - F_D(S)) \\
\end{align}
よって,dg(S)/dS = 0 から,
\begin{align}
c + (h - v) F_D(S) - (p + r) (1 - F_D(S)) &= 0 \\
F_D(S) &= \frac{p + r - c}{h + p + r - v}
\end{align}
になる.2 階微分は
\begin{align}
\frac{d^2g(S)}{dS^2} &= (h - v) f_D(S) + (p + r) f_D(S) \\
&= (h - v + p + r) f_D(S)
\end{align}
である.したがって,コスト関数 g(S) は凸関数であり,1 階微分が 0 になる点は最小値を与える.
コスト関数 g(S) を最小化するための最適発注量 S^* は
S^* = F_D^{-1}\left(\frac{p + r - c}{h + p + r - v}\right)
となる.ここで,F_D^{-1} は需要 D の累積分布関数の逆関数である.
新聞売り子問題において,より一般的に,在庫超過費用(overage cost)と 在庫不足費用(underage cost)を考慮する.
C_o = h + c - v, \quad C_u = p + r - c
臨界率
在庫超過費用 C_o は,在庫が余ったときのコストである.1 部の在庫超過に対し,保管コストと仕入れコストが発生するが,残存価額が得られないため,C_o = h + c - v となる.
在庫不足費用 C_u は,在庫が不足したときのコストである.1 部の在庫不足に対し,欠品コスト p と失われた販売機会の利益 r - c が発生するため,C_u = p + r - c となる.
従って,
\begin{align*}
S^* &= F_D^{-1}\left(\frac{p + r - c}{h + p + r - v}\right)\\
&= F_D^{-1}\left(\frac{C_u}{C_o + C_u}\right)
\end{align*}
が得られる.
初期在庫を考慮した新聞売り子問題*
新聞売り子の初期在庫を I とする.I \leq S^* の場合,最適発注量は S^* - I となる.すなわち,在庫量を S^* にすればよい.
また,g(S) は凸関数であり,発注量 S \geq 0 であるため,I > S^* の場合,何も発注しないことが最適である.
したがって,最適発注量は
Q =
\begin{cases}
S^* - I, & \text{if } I \le S^*, \\[6pt]
0, & \text{if } I > S^* .
\end{cases}
となる.
このような発注方式を Base Stock Policy (BSP) と呼ぶ.BSP は,各期間の在庫量を観測し,在庫量が S^* に引き上げられるように発注する方式である.新聞売り子問題において,BSP は最適な方策であると知られている.
需要が離散分布に従う場合*
需要 D が離散分布に従う場合は,
g(S) = h \sum_{d=0}^{S} (S - d) f_D(d) + p \sum_{d=S}^{\infty} (d - S) f_D(d)
となる.
ここで,有限差分(finite difference)を用いて g(S) の最小化を行うことができる.
S \in \{0, 1, 2, \ldots\} に対し,
\Delta g(S) = g(S + 1) - g(S)
とおくと,
\begin{align}
\Delta g(S) &= h \sum_{d=0}^{S + 1} (S + 1 - d) f_D(d) + p \sum_{d=S + 1}^{\infty} (d - S - 1) f_D(d) \\
&\quad - \left( h \sum_{d=0}^{S} (S - d) f_D(d) + p \sum_{d=S}^{\infty} (d - S) f_D(d) \right) \\
&= h \sum_{d=0}^{S} f_D(d) - p \sum_{d=S + 1}^{\infty} f_D(d) \\
&= h F_D(S) - p (1 - F_D(S))
\end{align}
g(S) は凸関数であるため,\Delta g(S) \geq 0 となる最小の S を求めればよい.すなわち,
F_D(S) \geq \frac{p}{h + p}
を満たす最小の S を求めればよい.
まとめ
発注量 S ,需要の観測値 d に対するコストは g(S, d) = h(S - d)^+ + p(d - S)^+ である.
コストの期待値 g(S) は凸関数であり,dg(S)/dS = 0 から F_D(S^*) = p/(h + p) が得られる.
したがって,最適発注量は S^* = F_D^{-1}\left(\frac{p}{h + p}\right) である.
p/(h + p) を臨界率と呼ぶ.F_D(S^*) は欠品が発生しない確率,すなわちサービスレベルであるため,サービスレベルを臨界率に一致させる発注量が最適である.
在庫不足費用 p が増加すると S^* は増加し,在庫超過費用 h が増加すると S^* は減少する.
*仕入れ価格や残存価額を含めた場合も,在庫超過費用 C_o = h + c - v ,在庫不足費用 C_u = p + r - c とおけば,S^* = F_D^{-1}\left(\frac{C_u}{C_o + C_u}\right) と同じ形になる.
文献案内
Arrow ほか (1951年 ) は,新聞売り子問題を初めて定式化した.
Scarf (1959年 ) の論文では,(s, S) 方策が発注費用を考慮した複数期間の新聞売り子問題において最適であることを示している.
実際には,需要の分布が不明であることが多い.Huber ほか (2019年 ) は,データ駆動新聞売り子問題(Data-Driven Newsvendor Problem)に関する研究が行われている.
Qin ほか (2011年 ) は,新聞売り子問題に関する研究をレビューした.
練習問題
練習 17.2 (Newsboy) ある新聞売り子は,新聞1部の在庫不足費用が p = 50 円,新聞1部の在庫超過費用が h = 25 円である.需要 D が正規分布 N(300, 60^2) に従うとき,最適な発注量を求めよ.
練習 17.3 (Down Jackets) あるファッション会社では,毎年10月から12月にかけてダウンジャケットを販売している.ダウンジャケットの生産コストは 6000 円であり,販売価格は 10000 円である.売れ残ったダウンジャケットは 20% の価格でセール販売される.この会社の OR 部門は,6月に生産計画を立てる必要があり,このときに発注量を決定しなければならない.今年度のダウンジャケットの需要 D が正規分布 N(900, 60^2) に従うと予測している.また,欠品した場合の機会損失は 1000 円である.このとき,最適な発注量を求めよ.
練習 17.4 (Bento Shop) ある弁当屋では,1 個 500 円で弁当を仕入れ,1 個 800 円で販売している.売れ残った弁当を廃棄する場合,1 個当たり 10 円の廃棄費用がかかる. 弁当の需要 D は,正規分布 N(50, 8^2) に従うとする.このとき,最適な弁当の発注量を求めよ.
解答 17.1 . 弁当 1 個当たりの在庫超過費用 h ,在庫不足費用 p は以下のように求められる.
h = 10 + 500 = 510, \quad p = 800 - 500 = 300
したがって,S^* は以下のように求められる.
S^* = F_D^{-1}\left(\frac{p}{h + p}\right) = F_D^{-1}\left(\frac{300}{510 + 300}\right) \approx F_D^{-1}\left(0.37\right)
標準正規分布表を調べると,\Phi(z) = 0.37 のとき,z は約 -0.33 である.したがって,S^* は
S^* = -0.33 \times \sigma + \mu = -0.33 \times 8 + 50 \approx 47.36
となる.
Python を用いて計算すると,以下のようになる.
from scipy.stats import norm
# 正規分布のパラメータ
mu = 50
sigma = 8
# 在庫超過費用と在庫不足費用
h = 510
p = 300
# 臨界率
critical_ratio = p / (h + p)
# 最適発注量
S_star = norm.ppf(critical_ratio, loc= mu, scale= sigma)
print (f"Optimal order quantity S*: { S_star:.4f} " )
Optimal order quantity S*: 47.3530
Arrow, Kenneth J, Theodore Harris, と Jacob Marschak. 1951年. 「Optimal Inventory Policy」 . Econometrica 19 (3): 250.
Huber, Jakob, Sebastian Müller, Moritz Fleischmann, と Heiner Stuckenschmidt. 2019年. 「A data-driven newsvendor problem: From data to decision」 . Eur. J. Oper. Res. 278 (3): 904–15.
Qin, Yan, Ruoxuan Wang, Asoo J Vakharia, Yuwen Chen, と Michelle M H Seref. 2011年. 「The newsvendor problem: Review and directions for future research」 . Eur. J. Oper. Res. 213 (2): 361–74.
Scarf, Herbert. 1959年. The optimality of (S , s) policies in the dynamic inventory problem .