在庫最適化の例題
mean = 100.
sigma = 10.
normal = norm(mean,sigma) #需要量
price = 100. #販売価格
purchase = 70. #仕入れ値
demand = normal.rvs(100)
lostsales = 0.
salvage = 10. #バイトに半値で売るとき 50 :-)
loyality = 0.5
Q = 100 #発注量
QLB = 50
QUB = 130
retail_profit, seven_profit = [], []
for Q in range(QLB,QUB):
profit_retail, profit_seven, lostsales_cost, salvage_value = 0,0,0,0
for d in list(demand):
purchase_cost = purchase*Q
profit = price*min(Q,d)
net_margin = profit-purchase_cost #粗利益（通常会計 -> 本社利益と小売店利益は同じ）
profit_retail += (net_margin + salvage*max(Q-d,0) -lostsales*max(d-Q,0))*(1.-loyality) #loyality = phi (=supplierへの利益分配率）)
profit_seven += (net_margin+ salvage*max(Q-d,0)- lostsales*max(d-Q,0) )*loyality
lostsales_cost += lostsales*max(d-Q,0) #機会損出費用
salvage_value += salvage*max(Q-d,0) #売れ残り価値
retail_profit.append( profit_retail)
seven_profit.append( profit_seven)
#print(profit_retail, profit_seven, lostsales_cost, salvage_cost)
#最適な発注量を求める => プロット
# min_cost = np.max
retail_profit2, seven_profit2 = [], []
for Q in range(QLB,QUB):
profit_retail, profit_seven, lostsales_cost, salvage_value = 0,0,0,0
for d in list(demand):
purchase_cost = purchase*min(Q,d) #売れた分だけ仕入れ値を払う（コンビニ会計）
profit = price*min(Q,d)
net_margin = profit-purchase_cost #粗利益（コンビニ会計だとnet_marginが大きくなる）
profit_retail += net_margin*(1.-loyality) -purchase*max(Q-d,0)+salvage*max(Q-d,0) -lostsales*max(d-Q,0)*(1.-loyality) #売れ残りに対する原価を店舗利益から引く（コンビニ会計）
profit_seven += net_margin*loyality- lostsales*max(d-Q,0)*loyality
lostsales_cost += lostsales*max(d-Q,0) #機会損出費用
salvage_value += salvage*max(Q-d,0) #売れ残り価値
retail_profit2.append( profit_retail)
seven_profit2.append( profit_seven)
#print(profit_retail, profit_seven, lostsales_cost, salvage_cost)
cost_df = pd.DataFrame({"発注量":list(range(QLB,QUB)),"小売利益":retail_profit, "本社利益":seven_profit,
"小売利益（コンビニ会計）":retail_profit2, "本社利益（コンビニ会計）":seven_profit2})
cost_df.head()
fig = px.scatter(cost_df, x="発注量",y=["小売利益", "本社利益","小売利益（コンビニ会計）", "本社利益（コンビニ会計）" ])
plotly.offline.plot(fig);
予測と在庫管理の融合
従来の需要予測の研究と在庫管理の研究は別々に行われてきた． 予測では，非定常な需要を仮定し，在庫管理では定常な分布を仮定する場合が多い． これらの2つの（おそらくサプライ・チェインに対して最も重要かつ最も多くの研究が行われてきた）モデルを構築することは，永年の研究者たちの夢であったが，いまだに実務的で納得するものは出ていない．
従来の研究の多くは，以下の２つに分類される．
- 綺麗な解析ができるような（実務的にはかなり特殊な）仮定の元での理論研究
- 実際の問題をなんとか解決するための（理論的な背景はほとんど考慮されていない）ヒューリスティクス
以下では，理論的にもある程度きちんとした枠組みで，かつ実際問題を解決可能な方法について考える．方針は以下の通り．
予測を行い誤差の分布を適当な分布で近似する．
リード時間分の予測需要の合計と分布に基づく安全在庫の和が変動型の基在庫レベルになる．
需要が負にならないようにスライドして，基在庫レベル最適化を行い，安全在庫量の最適化を行う．
安全在庫と需要予測に基づく確定在庫を合わせることによって，動的な基在庫レベルを導く．
Integrating Forecasting and Inventory Management
Traditional research on demand forecasting and inventory management have been conducted separately. Forecasting often assumes a non-stationary demand, while inventory management assumes a stationary distribution. Building these two (probably the most important and most studied) models for the supply chain has been a dream of researchers for a long time, but nothing practical and convincing has been produced yet.
Most conventional research can be divided into two categories.
- Theoretical research based on assumptions that allow for clean analysis (which are quite specific from a practical standpoint)
- Heuristics to somehow solve real problems (with little or no consideration of theoretical background).
In the following, we will consider methods that are both theoretically sound and capable of solving practical problems. The policy is as follows.
Make a forecast and approximate the distribution of the error with a suitable distribution.
The sum of the total demand forecast for the lead time and the safety stock based on the distribution is the base stock level for the variable type.
The base inventory level is optimized by sliding the demand so that it does not become negative, and the amount of safety stock is optimized.
The dynamic base stock level is derived by combining the safety stock and the fixed stock based on the demand forecast.
Loading sample demand data
Load the demand data with promotion data added. Choose customers and products, and generate a daily demand series. The data covers the period from the beginning of 2019 to the end of 2020, and the data after October 2020 is used for validation.
The forecasting method is based on automatic machine learning. This part can be replaced by deep learning or Bayesian inference.
promo_df = pd.read_csv(folder+"promo.csv", index_col=0)
demand_df = pd.read_csv(folder+"demand_with_promo_all.csv")
c = "仙台市"
p = "C"
print(c,p)
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= promo_df, agg_period="1d", forecast_periods = 10)
horizon="2020/10/01"
best, result_df = automl(df, horizon, 5)
fig, error, result = predict_using_automl(df, future, best[0])
#plotly.offline.plot(fig);
train_idx, valid_idx = prepare_index(df, horizon)
print(len(valid_idx))
error = result[:len(df)].Label - df.demand
error[valid_idx].hist(density=True);
error[valid_idx]
fig, best_dist, best_fit_name, best_fit_params = best_distribution(error[valid_idx])
#plotly.offline.plot(fig);
print(best_fit_name, best_fit_params)
d = best_dist.rvs((5,10000))
pd.Series(d.sum(axis=0)).hist(density=True);
n = norm(loc=best_dist.mean()*5,scale=best_dist.std()*np.sqrt(5))
d2 = best_dist.rvs(10000)
pd.Series(d2).hist(density=True);
Lmax = 50
MaxDemand = [0.]
for L in range(1,min(Lmax,len(data))): #リード時間
df_error = pd.Series(error[valid_idx])
data = df_error.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
fig, hist = best_histogram(data)
#plotly.offline.plot(fig2);
MaxDemand.append( hist.ppf(0.95) )
MaxDemand2 = []
max_ = MaxDemand[0]
for d in MaxDemand:
max_ = max(max_, d)
if d < max_:
MaxDemand2.append(max_)
else:
MaxDemand2.append(d)
pd.Series(MaxDemand2).plot();
L = 5#リード時間
df_error = pd.Series(error[valid_idx])
data = df_error.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
pd.Series(data).hist();
fig2, best_dist2, best_fit_name2, best_fit_params2 = best_distribution(data)
plotly.offline.plot(fig2);
print(best_fit_name2, best_fit_params2)
n_samples = 10
n_periods = 1000
capacity = 10.
convergence = 1e-5
b = 100
h = 1
LT = L
critical_ratio = b/(b+h)
S = best_dist2.ppf(critical_ratio)
lb = best_dist.ppf(0.01)
print(S, lb )
demand = best_dist.rvs((n_samples,n_periods)) - lb #需要を分布が0以上になるようにシフトさせる．
S = S - lb*LT #基在庫レベルはリード時間内の定常需要だけ大きくする．
print(S)
n_samples = 10
n_periods = 1000
capacity = 1000.
convergence = 1e-5
b = 100
h = 1
LT = L
critical_ratio = b/(b+h)
S = best_dist2.ppf(critical_ratio)
lb = best_dist.ppf(0.01)
print(S, lb )
demand = best_dist.rvs((n_samples,n_periods)) - lb #需要を分布が0以上になるようにシフトさせる．
S = S - lb*LT #基在庫レベルはリード時間内の定常需要だけ大きくする．
print(S)
print("t: S dS Cost")
for iter_ in range(100):
dC, cost, I = base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)
S = S - .1*dC
print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
if dC**2<=convergence:
break
pd.DataFrame(I[0]).plot(); #在庫の推移の可視化
safety = S + lb*LT
print(safety)
future = pd.Series(result[len(df):].Label)
data = future.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
data
data + safety
future = result.Label[valid_idx]
data = future.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
DS = data + safety #dynamic base stock level
dem = df.demand[valid_idx]
I = np.zeros( len(dem)+1 )
I[0] = DS[0] #エシェロン在庫ポジション
for t, d in enumerate(dem):
I[t+1] = I[t] - d
order = max(DS[t+1] -I[t+1],0)
I[t+1] = I[t+1] + order
#print(t,d,I[t+1])
if t>=len(DS)-2:
break
pd.Series(I[:-3]).plot();
dem.plot();
demand_df = pd.read_csv(folder+"demand_normal.csv")
c = "仙台市"
p = "A"
print(c,p)
agg_period ="1d"
df, future = make_forecast_df(demand_df, customer=c, product=p,promo_df= None, agg_period="1d", forecast_periods = 0)
df.demand.max()
horizon="2019/10/01"
best, result_df = automl(df, horizon, 5)
fig, error, result = predict_using_automl(df, future, best[0])
plotly.offline.plot(fig);
train_idx, valid_idx = prepare_index(df, horizon)
print(len(valid_idx))
error = result.Label - df.demand
error[valid_idx].hist(density=True);
fig, best_dist, best_fit_name, best_fit_params = best_distribution(error[valid_idx])
plotly.offline.plot(fig);
print(best_fit_name, best_fit_params)
L = 3 #リード時間
df_error = pd.Series(error[valid_idx])
data = df_error.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
fig2, best_dist2, best_fit_name2, best_fit_params2 = best_distribution(data)
plotly.offline.plot(fig2);
print(best_fit_name2, best_fit_params2)
n_samples = 10
n_periods = 1000
capacity = 1000.
convergence = 1e-5
b = 100
h = 1
LT = L
critical_ratio = b/(b+h)
S = best_dist2.ppf(critical_ratio)
lb = best_dist.ppf(0.01)
print(S, lb )
demand = best_dist.rvs((n_samples,n_periods)) - lb #需要を分布が0以上になるようにシフトさせる．
S = S - lb*LT #基在庫レベルはリード時間内の定常需要だけ大きくする．
print(S)
print("t: S dS Cost")
for iter_ in range(100):
dC, cost, I = base_stock_simulation(n_samples, n_periods, demand, capacity, LT, b, h, S)
S = S - 10.*dC
print(f"{iter_}: {S:.2f} {dC:.3f} {cost:.2f}")
if dC**2<=convergence:
break
pd.DataFrame(I[0]).plot(); #在庫の推移の可視化
safety = S + lb*LT
print(safety)
#検証データでシミュレーション
future = result.Label[valid_idx]
data = future.rolling(window=L).sum()[L-1:].values #L日前から直前までの需要の合計
DS = data + safety #dynamic base stock level
dem = df.demand[valid_idx]
I = np.zeros( len(dem)+1 )
I[0] = DS[0] #エシェロン在庫ポジション
for t, d in enumerate(dem):
I[t+1] = I[t] - d
order = max(DS[t+1] -I[t+1],0)
I[t+1] = I[t+1] + order
#print(t,d,I[t+1])
if t>=len(DS)-2:
break
pd.Series(I[:-3]).plot();
n = 7
G = SCMGraph()
for i in range(n):
G.add_node(i)
G.add_edges_from([(0, 2), (1, 2), (2,4), (3,4), (4,5), (4,6)])
z = np.full(len(G),1.65)
h = np.array([1,1,3,1,5,6,6])
mu = np.array([200,200,200,200,200,100,100])
sigma = np.array([14.1,14.1,14.1,14.1,14.1,10,10])
LTUB = np.array([0,0,0,0,0,3,1], int)
ProcTime = np.array([6,2,3,3,3,3,3], int) # 動的最適化の例題と合わせるため、生産時間から保証リード時間上限を減じてある
best_cost, best_sol, best_NRT, best_MaxLI, best_MinLT = tabu_search_for_SSA(G, ProcTime, LTUB, z, mu, sigma, h, max_iter = 10, TLLB =1, TLUB =3, seed = 1)
print("最良値", best_cost)
print("最良解", best_sol)
print("正味補充時間", best_NRT)
print("最大補充リード時間", best_MaxLI)
print("最小保証リード時間", best_MinLT)
pos = G.layout()
stage_df, bom_df = make_df_for_SSA(G, pos, ProcTime, LTUB, z, mu, sigma, h, best_NRT, best_MaxLI, best_MinLT)
stage_df.to_csv(folder_bom + "stage_ex1.csv")
bom_df.to_csv(folder_bom + "bom_ex1.csv")
fig, G, pos = draw_graph_for_SSA_from_df(stage_df, bom_df)
#fig.show()
nx.draw(G, pos=pos)
stage_df = pd.read_csv(folder_bom + "stage_ex1.csv")
bom_df = pd.read_csv(folder_bom + "bom_ex1.csv")
best_cost, stage_df, bom_df, fig = solve_SSA(stage_df, bom_df)
stage_df
cost, stage_df, I, cost_list, phi_list = periodic_inv_opt_fit_one_cycle(stage_df, bom_df, max_iter = 1000, n_samples = 10,
n_periods = 100, seed = 1, lr_find=False, max_lr =6.0, moms=(0.85,0.95))
