Lotsizing Optimzation using Gurobi & OptSeq

簡易システムの紹介ビデオ

ロットサイズ最適化システム OptLot

YouTubeVideo("0li4ma3EYs8")

多段階動的ロットサイズ決定モデル

Pochet--Wolsey 1991 Management Science Vol. 37 53--67 参照

ここでは,多段階にわたって製造を行うときのロットサイズ決定問題を考える.

標準定式化

集合

  • $\{1..T\}$: 期間の集合
  • $P$ : 品目の集合(完成品と部品や原材料を合わせたものを「品目」と定義する。)
  • $K$ : 生産を行うのに必要な資源(機械,生産ライン,工程などを表す)の集合
  • $P_k$ : 資源 $k$ で生産される品目の集合.
  • $Parent_p$ : 部品展開表における品目(部品または材料)$p$ の親品目の集合.言い換えれば,品目 $p$ から製造される品目の集合.

パラメータ

  • $T$: 計画期間数;期を表す添え字を $1,2,\cdots,t,\cdots,T$ と記す.
  • $f_t^p$ : 期 $t$ に品目 $p$ に対する段取り替え(生産準備)を行うときの費用(段取り費用)
  • $g_t^p$ : 期 $t$ に品目 $p$ に対する段取り替え(生産準備)を行うときの時間(段取り時間)
  • $c_t^p$ : 期 $t$ における品目 $p$ の生産変動費用
  • $h_t^p$ : 期 $t$ から期 $t+1$ に品目 $p$ を持ち越すときの単位あたりの在庫費用
  • $d_t^p$ : 期 $t$ における品目 $p$ の需要量
  • $\phi_{pq}$ : $q \in Parent_p$ のとき, 品目 $q$ を $1$ 単位製造するのに必要な品目 $p$ の数 ($p$-units); ここで, $p$-unitsとは,品目 $q$ の $1$単位と混同しないために導入された単位であり, 品目 $p$ の $1$単位を表す.$\phi_{pq}$ は,部品展開表を有向グラフ表現したときには,枝の重みを表す.
  • $M_t^k$ : 期 $t$ における資源 $k$ の使用可能な生産時間の上限. 定式化では,品目 $1$単位の生産時間を $1$単位時間になるようにスケーリングしてあるものと仮定しているが, プログラム内では単位生産量あたりの生産時間を定義している.
  • $UB_t^p$ : 期 $t$ における品目 $p$ の生産時間の上限. 品目 $p$ を生産する資源が $k$ のとき,資源の使用可能時間の上限 $M_t^k$ から段取り替え時間 $g_t^p$ を減じたものと定義される.

変数

  • $x_t^p$(x): 期 $t$ における品目 $p$ の生産量
  • $I_t^p$(inv) : 期 $t$ における品目 $p$ の在庫量
  • $y_t^p$(y): 期 $t$ に品目 $p$ に対する段取りを行うとき $1$, それ以外のとき $0$ を表す $0$-$1$ 変数

上の記号を用いると、多段階ロットサイズ決定モデルは,以下のように定式化できる.

$$ \begin{array}{ l l l } minimize & \sum_{t=1}^T \sum_{p \in P} \left( f_t^p y_t^p + c_t^p x_t^p + h_t^p I_t^p \right) & \\ s.t. & \ \ I_{t-1}^p +x_t^p = d_t^p+ \sum_{q \in Parent_p} \phi_{pq} x_t^q +I_t^p & \forall p \in P, t=1,\cdots,T \\ & \sum_{p \in P_k} x_t^p +\sum_{p \in P_k} g_t^p y_t^p \leq M_t^k & \forall k \in K, t=1,\cdots,T \\ & x_t^p \leq UB_t^p y_t^p & \forall p \in P, t=1,\cdots,T \\ & I_0^p =0 & \forall p \in P \\ & x_t^p,I_t^p \geq 0 & \forall p \in P, t=1,\cdots,T \\ & y_t \in \{0,1\} & \forall t=1,\cdots,T \end{array} $$

上の定式化で,最初の制約式は,各期および各品目に対する在庫の保存式を表す. より具体的には,品目 $p$ の期 $t-1$ からの在庫量 $I_{t-1}^p$ と生産量 $x_t^p$ を加えたものが, 期 $t$ における需要量 $d_t^p$,次期への在庫量 $I_t^p$, および他の品目を生産するときに必要な量 $\sum_{q \in Parent_p} \phi_{pq} x_t^q$ の和に等しいことを表す.

2番目の制約は, 各期の生産時間の上限制約を表す. 定式化ではすべての品目の生産時間は, $1$ 単位あたり$1$ 時間になるようにスケーリングしてあると仮定していたが,実際問題のモデル化の際には, 品目 p を $1$ 単位生産されるときに,資源 $r$ を使用する時間を用いた方が汎用性がある.

3番目の式は,段取り替えをしない期は生産できないことを表す.

エシェロン在庫を用いた定式化

品目間の親子関係だけでなく,先祖(部品展開表を表す有向グラフを辿って到達可能な品目の集合)を導入しておく.

集合

  • $Ancestor_{p}$: 品目 $p$ の先祖の集合. 親子関係を表す有向グラフを辿って到達可能な点に対応する品目から構成される集合. 品目 $p$ 自身は含まないものとする.

パラメータ

  • $\rho_{pq}$: $q \in Ancestor_p$ のとき,品目 $q$ を $1$ 単位生産するのに必要な製品 $p$ の量
  • $H_t^p$: 期 $t$ における品目 $p$ のエシェロン在庫費用; 品目 $p$ を生産することによって得られた付加価値に対する在庫費用を表す。品目 $p$ を製造するのに必要な品目の集合を $Child_p$ としたとき,以下のように定義される. $$ H_t^p =h_t^p -\sum_{q \in Child_p} h_t^q \phi_{qp} $$

変数

  • $E_{t}^p$: 期 $t$ における品目 $p$ のエシェロン在庫量;自分と自分の先祖の品目の在庫量を合わせたものであり, 以下のように定義される. $$ E_t^p = I_t^p +\sum_{q \in Ancestor_p} \rho_{pq} I_t^q $$

上の記号を用いると,エシェロン在庫を用いた多段階ロットサイズ決定モデルの定式化は,以下のようになる. $$ \begin{array}{ l l l } minimize & \sum_{t=1}^T \sum_{p \in P} \left( f_t^p y_t^p + c_t^p x_t^p + H_t^p E_t^p \right) & \\ s.t. & \ \ E_{t-1}^p +x_t^p - E_t^p = d_t^p + \sum_{q \in Ansestor_p} \rho_{pq} d_t^q & \forall p \in P, t=1,\cdots,T \\ & \sum_{p \in P_k} x_t^p +\sum_{p \in P_k} g_t^p y_t^p \leq M_t^k & \forall k \in K, t=1,\cdots,T \\ & E_t^p \geq \sum_{q \in Parent_p} \phi_{pq} E_t^q & \forall p \in P, t=1,\cdots,T \\ & x_t^p \leq UB_t^p y_t^p & \forall p \in P, t=1,\cdots,T \\ & E_0^p =0 & \forall p \in P \\ & x_t^p,E_t^p \geq 0 & \forall p \in P, t=1,\cdots,T \\ & y_t \in \{0,1\} & \forall t=1,\cdots,T \end{array} $$

上の定式化で,最初の制約は,各期および各品目に対するエシェロン在庫の保存式を表す. より具体的には,品目 $p$ の期 $t-1$ からのエシェロン在庫量 $E_{t-1}^p$ と生産量 $x_t^p$ を加えたものが, 期 $t$ における品目 $p$ の先祖 $q$ の需要量を品目 $p$ の必要量に換算したものの合計 $\sum_{q \in Ansestor_p} \rho_{pq} d_t^q$ と、 自身の需要量 $d_{t}^p$ と、次期へのエシェロン在庫量 $E_t^p$ の和に等しいことを表す. 2番目の制約は,各期の生産時間の上限制約を表す. 3番目の制約は,各品目のエシェロン在庫量が,その親集合の品目のエシェロン在庫量の合計以上であること(言い換えれば実需要量が負にならないこと)を規定する. 4番目の制約は,段取り替えをしない期は生産できないことを表す.

この定式化は1段階モデルと同じ構造をもつので、強化式を加えたり、施設配置定式化のような強い定式化に変形することが可能である。

以下のプログラムでは、実在庫量の上下限制約を付加している。エシェロン在庫モデルにおける3番目の制約のスラック変数が、実在庫量になっているので、 制約を以下のように書き換えた後で、実在庫量に対する上下限制約を付加すれば良い。

$$ E_t^p + I_t^p = \sum_{q \in Parent_p} \phi_{pq} E_t^q \ \ \ \forall p \in P, t=1,\cdots,T $$

施設配置定式化

ここでは,オリジナルと定式化と同じ $O(T)$ 個の $0$-$1$ 変数を用いて,より強い定式化を導く. 簡単のため1段階モデルを考える。品目を表す添え字 $p$ を付加し、在庫をレシェロン在庫に置き換えれば、多段階にも使うことができる。

我々の想定しているロットサイズ決定問題では品切れは許さないので, 期 $t$ の需要は,期 $t$ もしくはそれ以前の期で生産された量でまかなわれなければならない. そこで,期 $t$ の需要のうち,期 $s (s \leq t)$ で生産した量によってまかなわれた比率を表す変数 $X_{st}$ を導入する.

$X_{st}$ が $1$ のとき,期 $t$の需要量 $d_t$ が,期 $s$ で生産され,期 $t$ まで在庫される. 需要 $1$ 単位あたりの変動費用は,生産変動費用 $c_s$ と在庫費用 $\sum_{\ell=s}^{t-1} h_{\ell}$ の和になるので, $X_{st}$ の係数 $CH_{st}$ は,以下のように定義される. $$ CH_{st}= \left( c_s +\sum_{\ell=s}^{t-1} h_{\ell} \right) d_t $$

すべての需要は満たされなければならず, また $X_{st} >0$ のときには期 $s$ で生産しなければならないので, $y_s=1$ となる. ここで,$y_t$ はオリジナルの定式化における $0$-$1$ 変数で, 期 $t$ に生産をするときに $1$,それ以外のとき $0$ を表すことを思い起こされたい. よって,以下の定式化を得る.

$$ \begin{array}{ l l l } minimize & \sum_{s t: s \leq t} CH_{st} X_{st} + \sum_{t} f_t y_t & \\ s.t. & \ \ \sum_{s=1}^t X_{st}=1 & \forall t=1,\cdots,T \\ & X_{st} \leq y_s & \forall s \leq t, t=1,\cdots,T \\ & X_{st} \geq 0 & \forall s \leq t, t=1,\cdots,T \\ & y_t \in \{0,1\} & \forall t=1,\cdots,T \end{array} $$

この定式化は,期を表す点を需要地点(マーケット)ならびに施設(倉庫)の配置可能地点ととらえると, 施設配置問題に類似した問題になるので,施設配置定式化(facility location formulation)とよばれる.

オリジナルの定式化の変数との関係は以下のようになる。

各期の生産量: $$ x_s = \sum_{t=s}^T d_t X_{st} $$

在庫量: $$ I_t = \sum_{i=1}^t ( x_i -d_i ) $$

この関係を用いれば、オリジナルの定式化に含まれている付加制約を記述することができる。

他にも、生産をする際のロットサイズが固定されている場合(フル容量生産)に対しては、別の強化された定式化が可能である。 また、順序依存の段取りや複数の生産資源制約や生産ラインへの割り当てを考慮した定式化も可能であるが、問題の規模が大きくなると計算が困難になる。 実際問題の制約に応じて、カスタマイズされた定式化を行うことが望ましい。

# prod_df = pd.read_csv(folder + "Prod.csv", index_col="name")
# #prod_df.head()

# forecast_df = pd.read_csv(folder + "forecast_yhat.csv",index_col=0)
# #forecast_df.head()

# #予測を製品ごと(本来ならば工場ごと)に集約
# demand_ = forecast_df.groupby(by="prod").sum()
# demand = demand_.iloc[:,2:].values
# _, T = demand.shape 
# print("T=",T)

# #2段階モデルを仮定;各段階には1つの資源制約と仮定
# num_stages = 2
# maintenance_cycle = 100 #資源が使用不可になる周期;この数で割り切れると休止
# #num_resources = [2,3]
# #assert num_stages == len(num_resources) #段階数と資源数のリストの長さは同じであることの確認

# prod_time_bound = (1,1)
# setup_time_bound =(3600,7200)
# prod_cost_bound = (100,300)
# setup_cost_bound =(10000,20000)

# units_bound =(1,1)

# num_parents = 5  # 1つの子製品から生成される親製品の数;distribution型の生産工程を仮定
# products = list(prod_df.index)
# num_prod = len(products)

# #原材料リストを生成
# raw_material_name=""
# raw_materials = [] 
# for i,p in enumerate(prod_df.index):
#     raw_material_name += str(p)
#     if (i+1)%num_parents == 0 or i==len(prod_df)-1:
#         #print(i,p, raw_material_name)
#         raw_materials.append(raw_material_name)
#         raw_material_name =""
# #親子関係を定義
# parent = defaultdict(list)
# child = {}
# for r in raw_materials:
#     parent[r] = list(r)
#     for p in parent[r]:
#         child[p] = r

# #生産工程の容量を計算
# average_demand = demand.sum()/T
# cycle_time= prod_df["lot_size"].mean()/prod_df["average_demand"].mean()
# capacity = average_demand*prod_time_bound[1]+setup_time_bound[1]*num_prod/cycle_time

# print("capacity=",capacity)
# #資源データフレーム
# #期ごとに容量が変化するモデル
# #print(capacity)
# name_, period_, capacity_ = [], [], []
# for s in range(num_stages):
#     for t in range(T):
#         name_.append( f"Res{s}")
#         period_.append(t)
#         if (t+1) % maintenance_cycle == 0:
#             capacity_.append( 0 )
#         else:
#             if s ==0:
#                 capacity_.append( capacity/2 )
#             else:
#                 capacity_.append( capacity)
# resource_df = pd.DataFrame(data={"name": name_, "period": period_, "capacity": capacity_})


# #部品展開表のデータフレーム生成
# bom_df = pd.DataFrame(data={"child": [child[p] for p in products],
#                          "parent": products,
#                          "units": [random.randint(units_bound[0], units_bound[1]) for i in range(num_prod)]
#                         })

# #生産情報データフレーム生成
# items = products+raw_materials
# num_item = len(items)
# production_df = pd.DataFrame(data={"name": items,
#                              "ProdTime": [random.randint(prod_time_bound[0], prod_time_bound[1]) for i in range(num_item)],
#                              "SetupTime": [random.randint(setup_time_bound[0], setup_time_bound[1]) for i in range(num_item)],
#                              "ProdCost": [random.randint(prod_cost_bound[0], prod_cost_bound[1]) for i in range(num_item)],
#                              "SetupCost": [random.randint(setup_cost_bound[0], setup_cost_bound[1]) for i in range(num_item)]
#                              })
# production_df.set_index("name", inplace=True)
# production_df.to_csv(folder+"production.csv")
# bom_df.to_csv(folder+"bom.csv")
# resource_df.to_csv(folder+"resource.csv")
#resource_df.head()
#bom_df.head()
#production_df

ロットサイズ決定問題を解く関数

2段階モデルを仮定し、最初の段階(原材料生成工程)における資源名を Res0、次の段階(完成品生産工程)における資源名をRes1とする。

原材料の在庫は許さないものとする。

親製品は子製品1単位から生成される。

引数:

  • prod_df : 製品データフレーム
  • production_df : 生産情報データフレーム
  • bom_df : 部品展開表データフレーム
  • forecast_df : 予測データフレーム
  • resource_df : 資源データフレーム

返値:

  • model : ロットサイズ決定モデルのオブジェクト(最適解の情報も mode.__data に含む); statusも返すべきか?
  • T : 計画期間数

lotsizing[source]

lotsizing(prod_df, production_df, bom_df, forecast_df, resource_df)

ロットサイズ決定問題を解く関数

lotsizing関数の使用例

prod_df = pd.read_csv(folder + "Prod_with_inventory.csv",index_col="name")
production_df = pd.read_csv(folder + "production.csv",index_col="name")
bom_df = pd.read_csv(folder + "bom.csv")
forecast_df = pd.read_csv(folder + "forecast_yhat.csv")
resource_df = pd.read_csv(folder + "resource.csv")
if False: #pulpだと時間がかかり nbdev_test_nbs でタイムアウトするので、実験時のみTrueにする。
    model, T = lotsizing(prod_df, production_df, bom_df, forecast_df, resource_df)

最適化結果から図とデータフレームを生成する関数

引数:

  • model : ロットサイズ決定モデル
  • T : 計画期間
  • production_df : 生産情報データフレーム
  • bom_df : 部品展開表データフレーム
  • resource_df : 資源データフレーム

返値:

  • violated : 需要満足条件を逸脱した製品と期と逸脱量を保存したデータフレーム
  • production : 生産量を保管したデータフレーム
  • inventory : 在庫量を保管したデータフレーム
  • fig_inv : 在庫量の推移を表した図オブジェクト
  • fig_capacity : 容量制約を表した図オブジェクト

show_result_for_lotsizing[source]

show_result_for_lotsizing(model, T, production_df, bom_df, resource_df)

最適化結果から図とデータフレームを生成する関数

if False: # modelが存在するときのみTrueにする。
    violated, production, inventory, fig_inv, fig_capacity = show_result_for_lotsizing(model, T, production_df, bom_df, resource_df)
#fig_inv.show()
Image("../figure/fig_inv_for_lotsizing.png")
#fig_capacity.show()
Image("../figure/fig_capacity_for_lotsizing.png")