はじめに
今回は数理最適化の練習として、多期間にわたる在庫モデルの問題をPythonを使って解いてみたいと思います。
制約条件が複雑化したときの定式化の参考になれば幸いです。
問題概要
参考例題としてMathworksのHPで公開されていたMATLAB用の例題を一部※使います。
※本当は全部解きたかったんですが、変数設定と制約条件が200近くあり、これをうまくForループで処理することができなかったため。そのうちやります。
問題は以下のようになります。
経時的コスト変動が予測できるさまざまな原料を使用した、ある期間の配合肥料の生産スケジュールを立てたいです。あらかじめ、肥料への需要はわかっているものとします。目的は、需要を満たしながら利益を最大化することです。コストは、原材料の購入と経時的な肥料の保管にかかります。
また、制約条件は以下になります。
①2つの肥料製品（Balanced, HighN）のそれぞれの含有栄養素の重量比
②原材料の名前と栄養素の重量比※原料のSandには栄養素は無いですが、必要に応じてほか原料を希釈して栄養素の必要な重量比になるよう調整します。
③製品需要（精度の高い需要予測ができている想定。単位はトン）
④製品それぞれの１トン当たりの価格（＄）
- Balanced：４００
- HighN：５５０
⑤月ごとの原材料の価格（＄/t）（価格予測が精度良くできている想定）
⑥各種コスト・制約
- 保管コスト：10$/t・月
- 倉庫容量：1000t
- 月々の最大生産量：1200t
⑦その他
各月で満たされなかった需要は失われるものとします。つまり、期間内に注文を満たすことができない場合、超過分の注文は次の期間に繰り越されません。
Google OR-toolで解く
今回は1~3月分の生産計画を立ててみたいと思います。
③の製品需要を見ると1月、2月は最大生産量（1200t）以下の需要量なのでその月に必要量を生産すれば需給が間に合いますが、3月に入ると急に需要が増え、合わせて1500トンの需要が発生します。
これに対応するために需要が少ない月に多めに生産して、在庫として持っておく必要があります。
以下ではこれをGoogle OR-toolを使って解いて生産計画を立てて見ようと思います。
まず固定の変数を設定します。
変数設定
#固定の変数を設定
inventory_cost_rate =10
inventory_capacity = 1000
production_capacity = 1200
blend_price = 400
highn_price = 550
次に月ごとの情報をdfとして保存しておきます。
（本当はdfを使ってforループを回して変数設定、制約条件の追加をしたかったのですが、エラーが解消できなかったので今回は手打ちでやっています）
月ごとの需要・価格情報をデータフレーム化
import pandas as pd
import numpy as np
demand_Balanced = [750, 800, 900, 850, 700, 700, 700, 600, 600, 550, 550, 550]#1~12月のBalanced需要量
demand_HighN = [300, 310, 600, 400, 350, 300, 200, 200, 200, 200, 200, 200]#1~12月のHighN需要量
rawcost_MAP = [350, 360, 350, 350, 320, 320, 320, 320, 320, 310, 310, 340]#1~12月のMAP原料価格
rawcost_Potash = [610, 630, 630, 610, 600, 600, 600, 600, 600, 600, 600, 600]#1~12月のPotash原料価格
rawcost_AN = [300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300]#1~12月のAN原料価格
rawcost_AS = [135, 140, 135, 125, 125, 125, 125, 125, 125, 125, 125, 125]#1~12月のAS原料価格
rawcost_TSP = [250, 275, 275, 250, 250, 250, 250, 240, 240, 240, 240, 240]#1~12月のTSP原料価格
rawcost_Sand = [80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80, 80]#1~12月のSand原料価格
df = pd.DataFrame(np.arange(96).reshape(12,8), columns=["demand_Balanced", "demand_HighN", "rawcost_MAP", "rawcost_Potash", "rawcost_AN", "rawcost_AS", "rawcost_TSP", "rawcost_Sand"], index=["1月", "２月", "3月", "4月", "5月", "6月", "7月", "8月", "9月", "10月", "11月", "12月"])
df["demand_Balanced"]=demand_Balanced
df["demand_HighN"]=demand_HighN
df["rawcost_MAP"]=rawcost_MAP
df["rawcost_Potash"]=rawcost_Potash
df["rawcost_AN"]=rawcost_AN
df["rawcost_AS"]=rawcost_AS
df["rawcost_TSP"]=rawcost_TSP
df["rawcost_Sand"]=rawcost_Sand
df = df.astype(float)
df.head(12)
dfの中身は以下です。
最後に線形計画問題をGoogle or-toolsを使って解きます。
or-toolsで線形計画問題を解く
#or-toolsの線形計画法のライブラリをインポート
from ortools.linear_solver import pywraplp
#ソルバーの設定（バックエンドはSCIP）。
solver = pywraplp.Solver.CreateSolver('SCIP')
#変数の設定。非負制約も入れておく。
infinity = solver.infinity()
#1月分の変数設定
p1_1= solver.NumVar(0.0, infinity, '')#1月のBalancedの生産量
p2_1= solver.NumVar(0.0, infinity, '')#1月のHighNの生産量
s1_1 = solver.NumVar(0.0, df["demand_Balanced"][0], '')#1月のBalanced販売量。需要量を超えない制約あり。
s2_1 = solver.NumVar(0.0, df["demand_HighN"][0], '')#1月のHighN販売量。需要量を超えない制約あり。
u1_1 = solver.NumVar(0.0, infinity, '')#1月のMAP原材料の使用量
u2_1 = solver.NumVar(0.0, infinity, '')#1月のPotash原材料の使用量
u3_1 = solver.NumVar(0.0, infinity, '')#1月のAN原材料の使用量
u4_1 = solver.NumVar(0.0, infinity, '')#1月のAS原材料の使用量
u5_1 = solver.NumVar(0.0, infinity, '')#1月のTSP原材料の使用量
u6_1 = solver.NumVar(0.0, infinity, '')#1月のSand原材料の使用量
production_1 = solver.NumVar(0.0, production_capacity, '')#1月の生産量合計
inventory1_1=solver.NumVar(0.0, inventory_capacity, '')#1月の在庫量Balanced
inventory2_1=solver.NumVar(0.0, inventory_capacity, '')#1月の在庫量HighN
inventory_1=solver.NumVar(0.0, inventory_capacity, '')#1月の在庫量合計
#2月分の変数設定
p1_2= solver.NumVar(0.0, infinity, '')#2月のBalancedの生産量
p2_2= solver.NumVar(0.0, infinity, '')#2月のHighNの生産量
s1_2 = solver.NumVar(0.0, df["demand_Balanced"][1], '')#2月のBalanced販売量。需要量を超えない制約あり。
s2_2 = solver.NumVar(0.0, df["demand_HighN"][1], '')#2月のHighN販売量。需要量を超えない制約あり。
u1_2 = solver.NumVar(0.0, infinity, '')#2月のMAP原材料の使用量
u2_2 = solver.NumVar(0.0, infinity, '')#2月のPotash原材料の使用量
u3_2 = solver.NumVar(0.0, infinity, '')#2月のAN原材料の使用量
u4_2 = solver.NumVar(0.0, infinity, '')#2月のAS原材料の使用量
u5_2 = solver.NumVar(0.0, infinity, '')#2月のTSP原材料の使用量
u6_2 = solver.NumVar(0.0, infinity, '')#2月のSand原材料の使用量
production_2 = solver.NumVar(0.0, production_capacity, '')#2月の生産量合計
inventory1_2=solver.NumVar(0.0, inventory_capacity, '')#2月の在庫量Balanced
inventory2_2=solver.NumVar(0.0, inventory_capacity, '')#2月の在庫量HighN
inventory_2=solver.NumVar(0.0, inventory_capacity, '')#2月の在庫量合計
#3月分の変数設定
p1_3= solver.NumVar(0.0, infinity, '')#3月のBalancedの生産量
p2_3= solver.NumVar(0.0, infinity, '')#3月のHighNの生産量
s1_3 = solver.NumVar(0.0, df["demand_Balanced"][2], '')#3月のBalanced販売量。需要量を超えない制約あり。
s2_3 = solver.NumVar(0.0, df["demand_HighN"][2], '')#3月のHighN販売量。需要量を超えない制約あり。
u1_3 = solver.NumVar(0.0, infinity, '')#3月のMAP原材料の使用量
u2_3= solver.NumVar(0.0, infinity, '')#3月のPotash原材料の使用量
u3_3 = solver.NumVar(0.0, infinity, '')#3月のAN原材料の使用量
u4_3 = solver.NumVar(0.0, infinity, '')#3月のAS原材料の使用量
u5_3 = solver.NumVar(0.0, infinity, '')#3月のTSP原材料の使用量
u6_3 = solver.NumVar(0.0, infinity, '')#3月のSand原材料の使用量
production_3 = solver.NumVar(0.0, production_capacity, '')#3月の生産量合計
inventory1_3=solver.NumVar(0.0, inventory_capacity, '')#3月の在庫量Balanced
inventory2_3=solver.NumVar(0.0, inventory_capacity, '')#3月の在庫量HighN
inventory_3=solver.NumVar(0.0, inventory_capacity, '')#3月の在庫量合計
#制約条件の設定
#1月分の制約条件
solver.Add(u1_1*0.11 + u2_1*0 + u3_1*0.35 + u4_1*0.21 + u5_1*0 + u6_1*0 == p1_1*0.1 + p2_1*0.2)#マテバラ：N
solver.Add(u1_1*0.48 + u2_1*0 + u3_1*0 + u4_1*0 + u5_1*0.46 + u6_1*0 == p1_1*0.1 + p2_1*0.1)#マテバラ：P
solver.Add(u1_1*0 + u2_1*0.6 + u3_1*0 + u4_1*0 + u5_1*0 + u6_1*0 == p1_1*0.1 + p2_1*0.1)#マテバラ：K
solver.Add(u1_1 +u2_1 + u3_1 + u4_1 + u5_1 + u6_1 == p1_1 + p2_1)#マテバラ：生産量
solver.Add(p1_1 + p2_1 == production_1)#マテバラ：生産量
solver.Add(s1_1 <=p1_1)#販売量の制約
solver.Add(s2_1 <=p2_1)#販売量の制約
solver.Add(inventory_1 <= inventory_capacity)#各月の在庫量はinventory_capacity以下
solver.Add(p1_1 + p2_1 <=production_capacity)#各月の生産量の制約
solver.Add(inventory1_1 + inventory2_1== inventory_1)#品目別の在庫の合計
solver.Add(p1_1 - s1_1== inventory1_1)#品目別在庫
solver.Add(p2_1 - s2_1== inventory2_1)#品目別在庫
#2月分の制約条件
solver.Add(u1_2*0.11 + u2_2*0 + u3_2*0.35 + u4_2*0.21 + u5_2*0 + u6_2*0 == p1_2*0.1 + p2_2*0.2)#マテバラ：N
solver.Add(u1_2*0.48 + u2_2*0 + u3_2*0 + u4_2*0 + u5_2*0.46 + u6_2*0 == p1_2*0.1 + p2_2*0.1)#マテバラ：P
solver.Add(u1_2*0 + u2_2*0.6 + u3_2*0 + u4_2*0 + u5_2*0 + u6_2*0 == p1_2*0.1 + p2_2*0.1)#マテバラ：K
solver.Add(u1_2 +u2_2 + u3_2 + u4_2 + u5_2 + u6_2 == p1_2 + p2_2)#マテバラ：生産量
solver.Add(p1_2 + p2_2 == production_2)#マテバラ：生産量
solver.Add(s1_2 <=p1_2 + inventory1_1)#販売量の制約。最初の月以降は前月の在庫も追加。
solver.Add(s2_2 <=p2_2 + inventory2_1)#販売量の制約。最初の月以降は前月の在庫も追加。
solver.Add(p1_2 - s1_2 + inventory1_1 == inventory1_2)#2月の品目別在庫。前月の残量追加。
solver.Add(p2_2 - s2_2 + inventory2_1 == inventory2_2)#2月の品目別在庫。前月の残量追加。
solver.Add(inventory1_2 + inventory2_2== inventory_2)#2月在庫の合計
solver.Add(inventory_2 <= inventory_capacity)#各月の在庫量はinventory_capacity以下
solver.Add(p1_2 + p2_2 <=production_capacity)#各月の生産量の制約
#3月分の制約条件
solver.Add(u1_3*0.11 + u2_3*0 + u3_3*0.35 + u4_3*0.21 + u5_3*0 + u6_3*0 == p1_3*0.1 + p2_3*0.2)#マテバラ：N
solver.Add(u1_3*0.48 + u2_3*0 + u3_3*0 + u4_3*0 + u5_3*0.46 + u6_3*0 == p1_3*0.1 + p2_3*0.1)#マテバラ：P
solver.Add(u1_3*0 + u2_3*0.6 + u3_3*0 + u4_3*0 + u5_3*0 + u6_3*0 == p1_3*0.1 + p2_3*0.1)#マテバラ：K
solver.Add(u1_3 +u2_3 + u3_3 + u4_3 + u5_3 + u6_3 == p1_3 + p2_3)#マテバラ：生産量
solver.Add(p1_3 + p2_3 == production_3)#マテバラ：生産量
solver.Add(s1_3 <=p1_3 + p1_1 + p1_2 -s1_1 - s1_2)#販売量の制約
solver.Add(s2_3 <=p2_3 +p2_1 + p2_2 - s2_1 - s2_2)#販売量の制約
solver.Add(p1_3 - s1_3 + inventory1_2 == inventory1_3)#3月の品目別在庫。前月の残量追加。
solver.Add(p2_3 - s2_3 + inventory2_2 == inventory2_3)#3月の品目別在庫。前月の残量追加。
solver.Add(inventory1_3 + inventory2_3== inventory_3)#3月在庫の合計
solver.Add(inventory_3 <= inventory_capacity)#各月の在庫量はinventory_capacity以下
solver.Add(p1_3 + p2_3 <=production_capacity)#各月の生産量の制約
#各費用計算
inventory_cost = (inventory_1+inventory_2+inventory_3)*inventory_cost_rate
sales = blend_price*(s1_1+s1_2+s1_3) + highn_price*(s2_1 + s2_2 + s2_3)
raw_costs = u1_1*df["rawcost_MAP"][0] + u2_1*df["rawcost_Potash"][0] + u3_1*df["rawcost_AN"][0] + u4_1*df["rawcost_AS"][0] + u5_1*df["rawcost_TSP"][0] + u6_1*df["rawcost_Sand"][0]\
+u1_2*df["rawcost_MAP"][1] + u2_2*df["rawcost_Potash"][1] + u3_2*df["rawcost_AN"][1] + u4_2*df["rawcost_AS"][1] + u5_2*df["rawcost_TSP"][1] + u6_2*df["rawcost_Sand"][1]\
+u1_3*df["rawcost_MAP"][2] + u2_3*df["rawcost_Potash"][2] + u3_3*df["rawcost_AN"][2] + u4_3*df["rawcost_AS"][2] + u5_3*df["rawcost_TSP"][2] + u6_3*df["rawcost_Sand"][2]
#目的関数の設定。収益とコストからなる利益を最大化する。
solver.Maximize(sales - inventory_cost - raw_costs)
#ソルバーで設定した混合整数問題を解く。
solver.Solve()
結果確認
Objective().Value()で目的関数の値、solution_value()で目的関数が最大となるときの説明変数を呼び出すことができます。
例えば以下のようなコードを実行すると計算結果を確認できます。
結果確認（一部）
print('Solution:')
print('Objective value =', solver.Objective().Value())#722982.6562995873
print('p1_1 =', p1_1.solution_value())#847.8260869565217
3ヶ月分の各変数を出力してデータフレームにまとめた結果は以下になります。（ちょっと小さいですが、拡大すると各月の変数の結果が見れます）
最初に予測したとおり、3月の需要に対応するため、1月2月で作りだめをするように生産計画が立てられていることがわかりました。
また、3月分は1月の在庫を0からスタートすると、どう頑張ってもすべての需要を満たすことができないのですが、結果として利益最大化のために単価の大きなHighNの需要は100％満たすように作りだめが行われ、Balancedの需要はできる限り満たすが一部不足分が失われる結果になりました。
このように設定した制約条件や最大化したい目的関数の設定がきちんと反映された生産計画が立てられていることがわかりました。
以下で2つの製品Balanced, HighNについて各月の最適化結果を確認してみます。
・1~3月の生産量
・1~3月の販売量
・1~3月の在庫量
おわりに
今回は3ヶ月分の毎月変動する需要と原料価格をもとに利益が最大化するように各月の生産量・原料使用量・在庫量の使用量の線形計画問題を解きました。
実際にはもう少し細かく複雑な制約が現場にはありますが、複数月にまたがる生産計画を人間の頭で考えることは難しいので、定式化できるものについては数理最適化手法を使って効率化していきたいですね。
また、最近ではGoogle or-toolsのような使いやすい数理最適化ライブラリも出てきていますので数理最適化の実験を誰でも無料で簡単にできるようになってきています。
ある意味、現代は数理最適化を使う人間側の制約条件が緩和され、これまでよりも何かをする際に複数の選択肢を持てる、そんな時代になって来たように思えます。
世の中をより良くするために数理最適化がこれまで以上に活用されて行くことを願っています。