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

2023年3月5日2023年10月16日

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

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

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

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

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

1. なぜ行列演算を使うのか？
2. 今回使用する需要データ
3. 行列演算を使う箇所
4. シミュレーション結果をグラフで表示する
5. まとめ

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

今回使用する需要データ

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

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

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

行列演算を使う箇所

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

前提条件を設定する

まずは５つの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ごとの在庫補充目標量を二次元の行列に格納します。

これは７５×５の行列になります。

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

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

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

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

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

この行列は７５×（リードタイム－１）×SKU数の三次元行列になりますが、リードタイムはすべてのSKUの中で一番長い値とします。

ld_max = ld_list.max() # 一番長いリードタイムを求める
back_log_matrix = np.zeros((75, ld_max - 1, items)) # 日々のバックオーダー数、すべてゼロ、行列サイズは日数×（リードタイムー１）×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] # （リードタイムー１日）前までに発注したバックオーダー数
            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)) # 日々のバックオーダー数、すべてゼロ、行列サイズは日数×（リードタイムー１）×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] # （リードタイムー１日）前までに発注したバックオーダー数
            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,
　　　　……….

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

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

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

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()