企業や産業システムでは、原料・製品が供給から需要へと流れていく過程で それらを貯蔵する機能が必要である。この貯蔵容量と、それを運用管理する 最適の方法を決定する。
これらの合計を最小にし、利益を最大にする場合、 いつ、どれだけ、問屋に発注すれば良いのかを考える。
\(Q\)は販売店の在庫とし、一日に売れる個数を\(P\)個とする。 在庫量が基準発注点\(Q1\)をわったとき、\(a\)個の発注をする。 商品の入荷は\(lt\)日後とする。 ここで、売上個数\(P\)は毎日平均\(AL\)個のポアソン分布に従い、 入荷遅れ\(lt\)は平均\(ET\)、標準偏差\(ST\)の正規分布に従うとする。 さて、\(Q1\)の最適な値は何か。
設定を以下の通りで定める。
tend <- 1000 # 期間（日）
#tend1 <- 150 # Q1=10と固定した時にシミュレートした期間（日）
#P <- rpois(1,AL) # 売れた個数（個）
a <- 100 # 一回の入荷量（個）
c1 <- 20 # 一日一個当たりの在庫費用（円/個・日）
c2 <- 10000 # 一回の発注に伴う費用（円/回）
c3 <- 500 # 品切れ一個あたりに伴う損失（円/個）
SP <- 7000 # 売値
PC <- 5000 # 仕入れ単価
AL <- 5 # 一日の平均売上個数（個）
ET <- 6 # 納品するまでの平均日数（日）
ST <- 2 # 納品されるまでの日数の標準偏差
k <- 30 # jの取りうる値
lt <- -1 # 商品が届くまでの時間（日）：平均ET、標準偏差STの正規分布に従う
SC <- 0 # i日までの在庫量の総数（個）
nn <- 0 # 発注回数（回）
LS <- 0 # 売りそこなったものの総数（個）
YK <- 0 # 売上個数の総数
Q <- rep(NA,tend) # 現在庫量（個）
#SC <- # i日までの在庫量の総数（個）
#nn <- # 発注回数（回）
#LS <- # 売りそこなったものの総数（個）
#YK <- # 売上個数の総数
Q1 <- rep(NA,k) # 発注基準個数（個）
RR <- rep(NA,k) # 総売り上げ
Q[1] <- a+5 # 初期の在庫数（個）
#lt <- -1 # 発注していない状態
SC <- Q[1] #一日目の在庫数の総数
nn <- 0 #一日目の発注回数は0
LS <- 0 #一日目に売りそこなったものの数
YK <- 0 #一日目は売った個数ゼロ発注基準個数を\(Q1＝10\)個とした時の在庫変化のシミュレーションを与える。
#わかりやすくするため、150日間のシミュレーションを行う。
tend1 <- 150
# Q1をj=2,つまりQ1=10と固定する。
j <- 2
Q1[j] <- j*5 # 発注基準個数（個）
Q <- rep(NA, tend1) # 現在庫量（個）
Q[1] <- a+5 # 初期の在庫数（個）
#グラフ設定
xlim <- c(0,tend1)
ylim <- c(-20,Q1[j]+a)
# Q1=10の時のtend期間での在庫変化シミュレーション
for(i in 2:tend1){
#●売り出し
P <-rpois(1,AL) # 売れた個数（個）
Q[i] <- Q[i-1]-P # 在庫は（昨日の在庫量）-（今日の売り上げ個数）
# ここで需要が供給を上回ってしまった場合のプログラムが必要
if(Q[i-1] > 0){ # 在庫があるとき
if(Q[i]<0){ # 昨日の在庫数が正で今日の在庫数が負の時
LS <- LS - Q[i] ; # 売りそこなった数が増える
SC <- SC ; # 在庫は無いので総数は昨日のまま
YK <- YK + Q[i-1]; # 売った数は昨日の在庫分全て
Q[i] <- 0 # 在庫は0としておく
}else{ # 昨日の在庫数が正で今日の在庫数が正の時
LS <- LS ; # 売りそこなった数無し
SC <- SC + Q[i]; # 在庫の総数に今日の在庫数が加算される
YK <- YK + P # 売上総数に今日の売り上げ分が加算される
}
}else{ # 在庫が無い時
LS <- LS-Q[i] ; YK <- YK; # 売りそこなった数が増え、売上無し
if(Q[i]<0){ # 昨日の在庫数が負で、今日の在庫数も負の時
SC <- SC; # 在庫は無いので総数は昨日のまま
Q[i] <- 0 # 在庫は0としておく
}else{ # 昨日の在庫数が正で、今日の在庫数が正の時
SC <- SC + Q[i] # 在庫の総数に加算される
}
}
#●在庫の補充
lt <- lt-1 # 今日から見て在庫が届くまでの日数は昨日より1日引いたもの。
if(lt == 0){
Q[i] <- Q[i]+a ; nn <- nn +1
# もし今日が入荷日であれば在庫は今日の個数＋a。発注回数が1増える。
}else{Q[i] <- Q[i] ; nn <- nn
# 入荷日でなければ在庫は今日の個数のまま。発注回数そのまま。
}
if((Q[i] <= Q1[j]) && (lt < 0)) { # もし在庫が基準値より無く、注文していなければ
lt <- ceiling(abs(rnorm(1,ET,ST))) # 正規分布に従って〇日遅れで入荷されますよ。
}else{
if(lt < 0){
lt <- lt + 1 # 在庫が十分にあって、注文していなければlt=-1で保ちます。
}else{
lt <- lt # 在庫がない状態で注文済みの時はそのまま
# 在庫有りで注文済みという場合が仮にあったとしてもそのままとしておく
}
}
# 在庫Qの変化をプロットする
# plotのフォント指定
par(family="HiraKakuProN-W3")
plot(Q[1:i], xlim=xlim, ylim=ylim, xlab="日にちi(日)", ylab="在庫量Q(個)", type="l"); par(new=T)
# 発注基準個数Q1=10の線
par(family="HiraKakuProN-W3")
curve(0*x+Q1[j], from = 0, to = tend1, xlim=xlim, ylim=ylim, xlab="", ylab="", col="#1e90ff", type="l")
# Q1=10のときの利益RRを求める
SC1 <- SC * c1 # 在庫費用の累計
SC2 <- nn * c2 # 発注費用の累計
SC3 <- LS * c3 # 品切れ損失の累計
RR[j] <- (SP-PC) * YK -SC1-SC2-SC3
}# Q1=10と固定した時の在庫変化シミュレーション終わり↑\(tend1\)日での在庫量\(Q\)の変化が見られる。青の直線は発注基準個数\(Q1=10\)の線。この線を下回った時に注文をするという仕組みである。
# Q1をjと固定する。
# Q1=5,10,…,k*5 パターンを考える。
# グラフ設定
xlim <- c(0,tend)
ylim <- c(-20, k*5+a)
for(j in 1:k){
lt <- -1
# 以下先ほどの命令と同様------------------------------------------------------------------------
Q1[j] <- j*5 # 発注基準個数（個）
Q <- rep(NA, tend) # 現在庫量（個）
Q[1] <- a+5 # 初期の在庫数（個）
# Q1=j*5 の時のtend期間での在庫変化シミュレーション
for(i in 2:tend){
#●売り出し
P <-rpois(1,AL) # 売れた個数（個）
Q[i] <- Q[i-1]-P # 在庫は（昨日の在庫量）-（今日の売り上げ個数）
# 元々の在庫の在る・無しで「売り損ない数・在庫数・売った数」の総数が変わる
if(Q[i-1] > 0){
if(Q[i]<0){
LS <- LS-Q[i] ;
SC <- SC ;
YK <- YK+Q[i-1];
Q[i] <- 0
}else{
LS <- LS ;
SC <- SC+Q[i];
YK <- YK+P
}
}else{
LS <- LS-Q[i] ;
YK <- YK ;
if(Q[i]<0){
SC <- SC;
Q[i] <- 0
}else{SC<- SC+Q[i]
}
}
#●在庫の補充
lt <- lt-1 # 今日から見て在庫が届くまでの日数は昨日より1日引いたもの。
if(lt == 0){
Q[i] <- Q[i]+a ; nn <- nn +1
# もし今日が入荷日であれば在庫は今日の個数＋a。発注回数が1増える。
}else{Q[i] <- Q[i] ; nn <- nn
# 入荷日でなければ在庫は今日の個数のまま。発注回数そのまま。
}
if((Q[i] <= Q1[j]) && (lt < 0)) { # もし在庫が基準値より無く、注文していなければ
lt <- ceiling(abs(rnorm(1,ET,ST))) # 正規分布に従って〇日遅れで入荷されますよ。
}else{
if(lt < 0){
lt <- lt + 1 # 在庫が十分にあって、注文していなければlt=-1で保ちます。
}else{
lt <- lt # 在庫がない状態で注文済みの時はそのまま
# 在庫有りで注文済みという場合が仮にあったとしてもそのままとしておく
}
}
# Q1=j*5のときの利益RRを求める
SC1 <- SC * c1 # 在庫費用の累計
SC2 <- nn * c2 # 発注費用の累計
SC3 <- LS * c3 # 品切れ損失の累計
RR[j] <- (SP-PC) * YK -SC1-SC2-SC3
}
# Q1=j*5と固定した時の在庫変化シミュレーション終わり↑-------------------------------------------
# 在庫量変化の相違をグラフに表わす。
par(family="HiraKakuProN-W3")
plot(Q, xlim=xlim, ylim=ylim, xlab="日にちi(日)", ylab="在庫量Q(個)", type="l")
; par(new=T)
# 発注基準個数Q1の線
par(family="HiraKakuProN-W3")
curve(0*x+Q1[j], from = 0, to = tend, xlim=xlim, ylim=ylim, xlab="", ylab="", col="#1e90ff", type="l")
}発注基準個数\(Q1\)が上がるにつれて、在庫が足りなくなるということが無くなる。しかし、在庫がありすぎても在庫を管理する費用がかかる。
そこで、発注基準個数\(Q1\)と利益\(RR\)の関係を求めよう。
先ほどと同様にシミュレートし、得た値を表示する。
“Q1” : 発注基準個数（個） “YK” : 売上個数の累計（個） “LS” : 品切れのため、売り損なった個数の累計（個） “nn” : 発注回数（回） “SC1”: 在庫費用の累計 “SC2”: 発注費用の累計 “SC3”: 品切れ損失の累計 “RR” : 総売り上げ
を表す。
## [1] "Q1" "YK" "LS" "nn" "SC1" "SC2" "SC3" "RR"## [1] 5 3755 1194 37 720800 370000 597000 6001245
## [1] 10 4061 871 40 789040 400000 435500 6362020
## [1] 15 4050 889 40 810420 400000 444500 6712103
## [1] 20 4323 727 43 819260 430000 363500 7063754
## [1] 25 4566 448 45 908260 450000 224000 7401690
## [1] 30 4685 330 46 936540 460000 165000 7703509
## [1] 35 4851 295 48 991680 480000 147500 7937367
## [1] 40 4887 130 48 1063380 480000 65000 8092306
## [1] 45 4946 58 49 1181800 490000 29000 8149006
## [1] 50 5034 31 50 1299300 500000 15500 8145755
## [1] 55 5132 26 51 1308580 510000 13000 8123118
## [1] 60 4876 24 49 1466000 490000 12000 7988382
## [1] 65 5024 7 50 1527400 500000 3500 7899967
## [1] 70 4949 1 49 1661760 490000 500 7824481
## [1] 75 4907 0 49 1773600 490000 0 7710507
## [1] 80 5010 0 50 1824900 500000 0 7637355
## [1] 85 5139 0 51 1930540 510000 0 7514446
## [1] 90 5073 0 51 2078680 510000 0 7439540
## [1] 95 4914 0 49 2223720 490000 0 7315893
## [1] 100 4978 0 50 2237320 500000 0 7247692
## [1] 105 4938 0 50 2349800 500000 0 7138541
## [1] 110 4964 0 50 2450420 500000 0 7032785
## [1] 115 5078 0 51 2515080 510000 0 6932476
## [1] 120 5035 0 51 2690460 510000 0 6833238
## [1] 125 4970 0 50 2725260 500000 0 6736912
## [1] 130 4852 0 49 2906200 490000 0 6651883
## [1] 135 5049 0 51 2991900 510000 0 6539819
## [1] 140 5184 0 52 3028280 520000 0 6456680
## [1] 145 4950 0 50 3153060 500000 0 6344632
## [1] 150 4983 0 51 3226220 510000 0 6258201上のグラフを元に\(Q1\)と\(RR\)の関係を表すグラフを作成する。
ここではデータを3次曲線により近似する。 \(Q1\)を無限大に飛ばすと、利益は0に近づくが3次曲線の場合 無限大に行くので、本来この近似のやり方は好ましくない。 区間や領域を絞って回帰モデル化する局所回帰でやるとよい。
グラフが最大値をとる点が最大利益の点であるので、その時の\(Q1\)が適切な発注基準個数となる。
# nls関数にて3次式で近似する。
result <- nls(RR ~ (a * Q1^3 + b * Q1^2 + c * Q1 + d), start = c(a=0, b=0, c=0, d=0))
predict.c <- predict(result)
xlim <- c(0,k*5)
ylim <- c(min(RR,predict.c),max(RR,predict.c))
for(i in 1:k){
# 発注基準個数Q1と利益RRの点をプロットする
par(family="HiraKakuProN-W3")
plot(Q1, RR, xlim=xlim, ylim=ylim, type="p", xlab="発注基準個数 Q1", ylab="利益 RR") ; par(new=T)
# 近似曲線を描く
par(family="HiraKakuProN-W3")
plot(Q1[1:i], predict.c[1:i], type="l", xlim=xlim, ylim=ylim, xlab="",ylab="")}利益が最大の時は以下の通り。
max(predict.c)## [1] 8065339グラフから適切な発注基準個数を求めることができる。
## elapsed
## 19.031