行列演算を使った適正在庫シミュレーションをPythonで実装する【複数SKU用】

2023年10月16日

前回の記事Pythonで適正在庫シミュレーションソフトを作ってみた【単一SKU用】では単一SKUについて適正在庫シミュレーションするソフトをPythonで作りました。

適正在庫はSKUごとに設定しないと意味がないので、このシミュレーションをすべてのSKUについて実行することになります。

しかし、何千、何万というSKUを取り扱っている会社にとっては、そのようなことは気の遠くなるような話しです。

 

そこで今回は、複数のSKUについて一挙にシミュレーションするソフトを作ってみます。

そのために行列を使います。

 

◆仕事や勉強の息抜きに。。。

なぜ行列演算を使うのか?

適正在庫を決めるにはリードタイム発注サイクル等を設定しておく必要があります。

単一SKUのシミュレーションであればリードタイム=3日、発注サイクル=2日のように決めればよいので、変数はスカラーでOKです。

スカラーとは1とか2とかいう普通の数字のことです。

これらのスカラーを入れる変数の箱として、例えばリードタイムであればld=3とか、発注サイクルであればoc=2とかいうように定義しておけばよかったわけです。

 

ところが複数SKUの場合、それぞれのSKUでリードタイムや発注サイクルが異なるかもしれません。

そのため、2つのSKUのシミュレーションを行うのであれば、ld=[3,5]とかoc=[2,7]とかいうようにベクトルを使います。

ベクトルは行数または列数が1、つまり一次元の行列のことです。

 

また単一SKUのシミュレーションでは、日々の在庫数を格納するのにベクトルを使いました。

75日間のシミュレーションをするのであれば、1×75のベクトルに日々の在庫数を格納したわけです。

これが複数SKUのシミュレーションになると二次元の行列に格納します。

2つのSKUについてシミュレーションする場合には、2×75の行列に日々の在庫数を格納します。

 

更に、発注残(バックオーダー数)を入れる箱は三次元の行列になります。

単一SKUで納品リードタイムが3日のシミュレーションをする場合、発注残は1日目の分と2日目の分があるため、75日間のシミュレーションをするのであれば2×75の行列に日々の発注残を格納しました。

これが5つのSKUについてシミュレーションするのであれば、2×75×5の三次元行列になるわけです。

 

今回使用する需要データ

今回はSKU数が5つの場合のシミュレーションを行います。

但し、SKU数が数千、数万になっても大丈夫なように拡張性も持たせるようにします。

このデータが’raw_data.csv’のファイルに入っています。

 

行列演算を使う箇所

需要データファイルを読み込んで学習データとテストデータに分けるところまでは単一SKUの場合と同じですので、それ以降を下記のように変更していきます。

 

前提条件を設定する

まずは5つのSKUそれぞれについて平均標準偏差を計算してベクトルに格納します。

av_list = [] # 各SKUの需要データの平均を入れる空のベクトル
sd_list = [] # 各SKUの需要データの標準偏差を入れる空のベクトル
for i in range(1, items+1):
  av = learn_data[:, i].mean() #平均の計算
  sd = learn_data[:, i].std(ddof = 1) #不偏標準偏差の計算
  av_list.append(av)
  sd_list.append(sd)

ここでitemsという変数を使っていますが、この変数には下記のコーディングによってシミュレーションするSKU数が入っています。

learn_df = raw_df.iloc[:45] # 前半45日分は学習データ
test_df = raw_df.iloc[45:] # 45日目以降はテストデータ
learn_data = learn_df.to_numpy()
test_data = test_df.to_numpy()
rows, columns = test_data.shape # テストデータの行数と列数を取得
days = rows # 行数はシミュレーション日数
items = columns - 1 # 列数-1はSKU数(最初の列は日付であるため)
days, items

次に各SKUについてリードタイム、発注サイクル、許容欠品率を設定します。

今回は簡単のため、すべてのSKUについて同じ値とします。

# 初期設定
ld_list = np.full(items, 3) # すべてのSKUのリードタイム=3日
oc_list = np.full(items, 2) # すべてのSKUの発注サイクル=2日
so_list = np.full(items, 0.05) # すべてのSKUの許容欠品率=5%

 

在庫補充目標量を計算する

次に安全在庫需要予測在庫在庫補充目標量をSKU別に計算して、ベクトルに格納しておきます。

# 安全在庫の計算
safety_stock_by_item = []
for i in range(1, items+1):
    safety_stock = sd_list[i-1] * np.sqrt(ld_list[i-1] + oc_list[i-1]) * norm.ppf(1 - so_list[i-1]) # 安全在庫の計算
    safety_stock_by_item.append(safety_stock)

# 需要予測在庫の計算
cycle_stock_by_item = []
for i in range(1, items+1):
    cycle_stock = np.sum(ld_list[i-1] + oc_list[i-1]) * av_list[i-1] # 需要予測在庫の計算
    cycle_stock_by_item.append(cycle_stock)

# 在庫補充目標量の計算
target_stock_by_item = []
for i in range(1, items+1):
    target_stock = safety_stock_by_item[i-1] + cycle_stock_by_item[i-1] # 在庫補充目標数の計算
    target_stock_by_item.append(target_stock)

 また、シミュレーションにおいては最初の在庫数がゼロだと絶対に欠品になってしまいますので、ゼロ以外の初期在庫数を設定しておきます。

ここでは在庫補充目標量を初期在庫とします。

# 初期在庫設定(在庫補充目標量)
ini_stock = target_stock_by_item
ini_stock

 在庫補充目標量の計算式については、こちらの記事で解説しています。

適正在庫を維持するための発注数の決め方をわかりやすく【定期発注方式の場合】

 

各リストを初期化する

次にシミュレーション結果を格納する各ベクトルまたは行列の要素を初期化しておきます。

まずは、日々の各SKUごとの在庫補充目標量を二次元の行列に格納します。

これは75×5の行列になります。

target_stock_list = []
for i in range(75):
    target_stock_list.append(target_stock_by_item)

次に日々の在庫数発注数入荷数を格納する行列の要素をすべてゼロに初期化します。

これらの行列はすべて75×5の二次元行列です。

current_stock_list = np.zeros((75, items)) # 日々の在庫数、すべてゼロ
order_qty_list = np.zeros((75, items)) # 日々の発注数、すべてゼロ
receive_qty_list = np.zeros((75, items)) # 日々の入荷数、すべてゼロ

次に日々の発注残(バックオーダー数)を格納する行列の要素をすべてゼロに初期化します。

この行列は75×(リードタイム-1)×SKU数の三次元行列になりますが、リードタイムはすべてのSKUの中で一番長い値とします。

ld_max = ld_list.max() # 一番長いリードタイムを求める
back_log_matrix = np.zeros((75, ld_max - 1, items)) # 日々のバックオーダー数、すべてゼロ、行列サイズは日数×(リードタイムー1)×SKU数

 

倉庫内在庫量を計算する

以上で初期設定は終了したので、日々の在庫数の推移を計算していきます。

単一SKUの場合と異なるのは、items(SKU数)でforループを回すところだけです。

for k in range(items):
    for i in range(75):
        if i == 0:
            current_stock_list[i, k] = ini_stock[k] - test_data[i , k + 1] + receive_qty_list[i, k] # 初日のみ初期在庫から当日の在庫数を計算
        else:
            current_stock_list[i, k] = current_stock_list[i - 1, k] - test_data[i , k + 1] + receive_qty_list[i, k] # 前日の在庫数から当日の在庫数を計算

 

適正発注量を計算する

これもitems(SKU数)でforループを回すところだけが、単一SKUの場合と異なります。

for k in range(items):
    for i in range(75):
    if i % oc_list[k] == 0:
            order_qty_list[i, k] = target_stock_list[i, k] - current_stock_list[i, k] - np.sum(back_log_matrix[i, :, k]) # 発注日のみ発注数を計算

 

バックオーダー数と入荷数を計算する

これもitems(SKU数)でforループを回すところだけが、単一SKUの場合と異なります。

for k in range(items):
    for i in range(75):
        for j in range(ld_max):
            if j < ld_max - 1:
                back_log_matrix[i, j, k] = order_qty_list[i-(j+1), k] # (リードタイムー1日)前までに発注したバックオーダー数
            else:
                receive_qty_list[i, k] = order_qty_list[i-(j+1), k] # リードタイム後に納入される入荷数

 

リストに格納された日々の在庫数推移を表示する

以上をまとめると次のようになります。

import numpy as np
import pandas as pd
import scipy.stats as st
from scipy.stats import norm
from matplotlib import pyplot as plt

raw_df = pd.read_csv('raw_data.csv')
raw_df['Date'] = pd.to_datetime(raw_df['Date'])

learn_df = raw_df.iloc[:45] # 前半45日分は学習データ
test_df = raw_df.iloc[45:] # 45日目以降はテストデータ
learn_data = learn_df.to_numpy()
test_data = test_df.to_numpy()
rows, columns = test_data.shape # テストデータの行数と列数を取得
days = rows # 行数はシミュレーション日数
items = columns - 1 # 列数-1はSKU数(最初の列は日付であるため)

av_list = [] # 各SKUの需要データの平均を入れる空のベクトル
sd_list = [] # 各SKUの需要データの標準偏差を入れる空のベクトル
for i in range(1, items+1):
    av = learn_data[:, i].mean() #平均の計算
    sd = learn_data[:, i].std(ddof = 1) #不偏標準偏差の計算
    av_list.append(av)
    sd_list.append(sd)

# 初期設定
ld_list = np.full(items, 3) # すべてのSKUのリードタイム=3日
oc_list = np.full(items, 2) # すべてのSKUの発注サイクル=2日
so_list = np.full(items, 0.05) # すべてのSKUの許容欠品率=5%

# 安全在庫の計算
safety_stock_by_item = []
for i in range(1, items+1):
    safety_stock = sd_list[i-1] * np.sqrt(ld_list[i-1] + oc_list[i-1]) * norm.ppf(1 - so_list[i-1]) # 安全在庫の計算
    safety_stock_by_item.append(safety_stock)

# 需要予測在庫の計算
cycle_stock_by_item = []
for i in range(1, items+1):
    cycle_stock = np.sum(ld_list[i-1] + oc_list[i-1]) * av_list[i-1] # 需要予測在庫の計算
    cycle_stock_by_item.append(cycle_stock)

# 在庫補充目標量の計算
target_stock_by_item = []
for i in range(1, items+1):
    target_stock = safety_stock_by_item[i-1] + cycle_stock_by_item[i-1] # 在庫補充目標数の計算
    target_stock_by_item.append(target_stock)

# 初期在庫設定(在庫補充目標量)
ini_stock = target_stock_by_item
ini_stock

# 各リストの初期化
target_stock_list = []
for i in range(75):
    target_stock_list.append(target_stock_by_item)
target_stock_list
target_stock_list = np.array(target_stock_list)
current_stock_list = np.zeros((75, items)) # 日々の在庫数、すべてゼロ
order_qty_list = np.zeros((75, items)) # 日々の発注数、すべてゼロ
receive_qty_list = np.zeros((75, items)) # 日々の入荷数、すべてゼロ
ld_max = ld_list.max() # 一番長いリードタイムを求める
back_log_matrix = np.zeros((75, ld_max - 1, items)) # 日々のバックオーダー数、すべてゼロ、行列サイズは日数×(リードタイムー1)×SKU数

for k in range(items):
    for i in range(75):
        for j in range(ld_max):
            if j < ld_max - 1:
                back_log_matrix[i, j, k] = order_qty_list[i-(j+1), k] # (リードタイムー1日)前までに発注したバックオーダー数
            else:
                receive_qty_list[i, k] = order_qty_list[i-(j+1), k] # リードタイム後に納入される入荷数

        if i == 0:
            current_stock_list[i, k] = ini_stock[k] - test_data[i , k + 1] + receive_qty_list[i, k] # 初日のみ初期在庫から当日の在庫数を計算
        else:
            current_stock_list[i, k] = current_stock_list[i - 1, k] - test_data[i , k + 1] + receive_qty_list[i, k] # 前日の在庫数から当日の在庫数を計算

        if i % oc_list[k] == 0:
            order_qty_list[i, k] = target_stock_list[i, k] - current_stock_list[i, k] - np.sum(back_log_matrix[i, :, k]) # 発注日のみ発注数を計算

各SKUについて日々の在庫推移が格納された行列current_stock_listを見てみると次のようになりました。

current_stock_list

array([[ 7.08612859e+03,  1.13929747e+03,  2.59929034e+04,
         1.45256316e+03,  1.20085662e+05],
       [ 6.92612859e+03,  1.11029747e+03,  1.98869034e+04,
         9.01563160e+02,  8.94316621e+04],
       [ 6.57112859e+03,  1.09429747e+03,  1.57889034e+04,
         5.72563160e+02,  5.61276621e+04],
       [ 5.56812859e+03,  1.14329747e+03,  1.86799034e+04,
         8.06563160e+02,  9.03046621e+04],
       [ 3.87112859e+03,  1.08829747e+03,  1.79659034e+04,
         7.88563160e+02,  8.39206621e+04],
       [ 3.04012859e+03,  1.10929747e+03,  2.52919034e+04,
    ……….

これでは直感的にわかりにくいので、グラフにしてみましょう。

 

シミュレーション結果をグラフで表示する

まずは5つのSKUについて、需要数と在庫数の推移を1つのグラフで描くと次のようになります。

x = test_data[:, 0]
plt.figure(figsize=(10, 6))
plt.plot(x, current_stock_list, label = raw_df.keys()[1:])
plt.plot(x, test_data[:, 1], label = 'Demand')
plt.legend()
plt.show

 

 

これでは各SKUの需要量のスケールが違い過ぎて、需要量の小さいSKUの傾向がよくわかりませんね。

そこで、SKUごとにグラフを分けて表示してみましょう。

x = test_data[:, 0]
fig = plt.figure(figsize = (10,20))

for i in range(1, items+1):
    ax = str(i)
    ax = fig.add_subplot(5, 1, i)
    ax.plot(x, current_stock_list[:, i-1], color = 'darkblue', label = 'Inventory')
    ax.plot(x, test_data[:, i], color = 'darkorange', label = 'Demand')
    ax.set_ylabel(raw_df.keys()[i])
    ax.legend(loc = 'upper left')

plt.show()

 

まとめ

多SKUの適正在庫シミュレーションはExcel VBAでも実装することができます。

【VBAサンプル付き】適正在庫シミュレーションをエクセルマクロで自動化(前編)

しかしSKUごとに結果をExcelシートに書き込んでいくため時間がかかります。

SKUごとにシートを分ける場合は、シートの数にも制限があります。

グラフも作るとなると尚更です。

 

その点、PythonはSKU数が増えても行列で一気に計算するので、非常に短時間でシミュレーションできます。

またプログラミングも簡単で、Excel VBAよりもいろいろな処理が柔軟にできると思いました。

今後、更に大きなデータでいろいろな分析をしていきます。

【SCM分析1】曜日/物流センター別の出荷傾向をグラフで可視化する