ある工場では平日に100名以上,土曜日に94名以上,日曜日に75名以上の作業員を必要としています.各作業員は5日間連続して出勤した後,2日間の休日を取得することになっています.シフトパターンは7種類あり,月曜日から5日間出勤するパターンをパターン1,火曜日からを2,...,日曜日からをパターン7とします.一度決定したシフトパターンはある期間は継続して利用し,しばらくは変更しないものとします.このとき,各曜日に出勤する作業員の合計が最低人数を満たし,工場全体で雇用しなけらばならない作業員の合計人数を最小にするようなシフト計画を立案せよ.
\( \) この問題を定式化してみよう.曜日を \(i~(1=\)月曜日\(,2=\)火曜日\(, \cdots,7=\)日曜日\()\) で表す.曜日 \( i \) の作業員必要人数を \( d_i \) とし,そのベクトルを \( {\bf{d}} = (d_1, d_2, \cdots, d_7)^{\rm T} \) とします.曜日 \( i \) のシフトパターン \( j~(j=1,2,\cdots,7) \) を \( a_{ij} \) として,\( a_{ij}=1 \) のときには曜日 \( i \) にシフトパターン \( j \) の作業員が出勤し,\( a_{ij}=0 \) のときには出勤しないとします.
このとき,シフトパターンは,\( i \) 行 \( j \) 列目の要素が \( a_{ij} \) となるような次の行列で表現できます.
上の行列 \( A \) の左端の列はシフトパターン1で勤務するの作業員の出勤形態です.一方で1行目の1になっている部分が月曜日に出勤しているパターンの一覧を意味していることに注意します.
また,シフトパターン \( j \) で出勤(つまり雇用)する作業員の数を \( x_j \) とし,そのベクトルを \( {\bf{x}} = (x_1, x_2, \cdots, x_7)^{\rm T} \) とします.このとき次の通りシフト計画問題を定式化することができます.
なお,\( {\bf{Z}}_{+} \) は 0 以上の整数を意味します.また,\( A {\bf{x}} \geq {\bf{d}} \) は次のように書くこともできます.
さらに行列形式で表記された上の制約式を分解して次のような連立不等式に書くことも可能です.
線形最適化問題と同様に,PuLP や SciPy を使えば,このようなシフト計画問題を Python で簡単かつ高速に解くことができます.
PuLP は Anaconda をインストールしただけではインストールされていないので,まずは PuLP をインストールします.PuLP のインストールは「Anaconda Prompt」または「Anaconda Powershell Prompt」で pip install pulp
を実行すれば良いでしょう.なお,学内の情報処理実習室では pip install pulp --user
のように --user
オプションを指定する必要があるかもしれません.また,インストールされているモジュールの一覧は pip list
で確認できます.Jupyter Lab や Jupyter Notebook のマジックコマンドを利用しても構いません.
まずは一つずつ変数や制約式を定義して PuLP でシフト計画問題を解いてみよう.
PuLP モジュールをインポートします.もしもエラーが表示されたら,準備を参考に,PuLP をインストールしよう.また,行列を利用するので,NumPyもインポートします.
モジュールのインポート
import pulp
import numpy as np
曜日ごとに必要な人数を格納するためのベクトルを定義します.
曜日ごとの必要人数の設定
d = np.array([100, 100, 100, 100, 100, 94, 75])
シフトパターンの行列を定義します.
シフトパターンの行列
A = np.array([
[1, 0, 0, 1, 1, 1, 1],
[1, 1, 0, 0, 1, 1, 1],
[1, 1, 1, 0, 0, 1, 1],
[1, 1, 1, 1, 0, 0, 1],
[1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1],
])
次に問題オブジェクトを生成します.問題オブジェクトは pulp.LpProblem()
関数で生成できます.このとき sense 引数には,最大化問題を解きたい場合には sense=pulp.LpMaximize
を,最小化問題を解きたい場合には sense=pulp.LpMinimize
を指定します.なお,sense 引数を省略した場合は,最小化問題として生成されます.
問題オブジェクトを作る
problem = pulp.LpProblem(name='シフト計画', sense=pulp.LpMinimize)
続いて変数を定義します.7個あるので面倒ですが,まずは一つひとつ定義します(後のステップではまとめて定義する方法も示します).変数は 0 以上の整数型であるので,lowBound=0
と cat=pulp.LpInteger
を指定していることに注意してください.
変数を定義する
x1 = pulp.LpVariable('x1', lowBound=0, cat=pulp.LpInteger)
x2 = pulp.LpVariable('x2', lowBound=0, cat=pulp.LpInteger)
x3 = pulp.LpVariable('x3', lowBound=0, cat=pulp.LpInteger)
x4 = pulp.LpVariable('x4', lowBound=0, cat=pulp.LpInteger)
x5 = pulp.LpVariable('x5', lowBound=0, cat=pulp.LpInteger)
x6 = pulp.LpVariable('x6', lowBound=0, cat=pulp.LpInteger)
x7 = pulp.LpVariable('x7', lowBound=0, cat=pulp.LpInteger)
次は目的関数です.7 個の変数の和になるので,ここではその通り指定します(これも後のステップで別の指定方法を示します).
目的関数を設定する
problem += x1 + x2 + x3 + x4 + x5 + x6 + x7
制約条件式を1行ごとに設定します.
制約条件を設定する
for i in range(len(d)):
problem += A[i,0]*x1 + A[i,1]*x2 + A[i,2]*x3 + A[i,3]*x4 + A[i,4]*x5 + A[i,5]*x6 + A[i,6]*x7 >= d[i]
上で定義した式と同じように問題を設定できているか確認します.
問題を表示する
print(problem)
シフト計画: MINIMIZE 1*x1 + 1*x2 + 1*x3 + 1*x4 + 1*x5 + 1*x6 + 1*x7 + 0 SUBJECT TO _C1: x1 + x4 + x5 + x6 + x7 >= 100 _C2: x1 + x2 + x5 + x6 + x7 >= 100 _C3: x1 + x2 + x3 + x6 + x7 >= 100 _C4: x1 + x2 + x3 + x4 + x7 >= 100 _C5: x1 + x2 + x3 + x4 + x5 >= 100 _C6: x2 + x3 + x4 + x5 + x6 >= 94 _C7: x3 + x4 + x5 + x6 + x7 >= 75 VARIABLES 0 <= x1 Integer 0 <= x2 Integer 0 <= x3 Integer 0 <= x4 Integer 0 <= x5 Integer 0 <= x6 Integer 0 <= x7 Integer
solve()
関数で問題を解きます.このとき計算のログが多く出力されます.オプションを付与するとログを出力しないように変更できます.
問題を解く
problem.solve()
# problem.solve(pulp.PULP_CBC_CMD(msg=0)) # ログを出力しない
1
問題を解けたかどうか念のために確認します.Optimal
と表示されたので正しく解けているようです.
結果の確認
print("Status : ", pulp.LpStatus[problem.status])
Status : Optimal
最後に結果を表示します.月曜日から出勤するシフトパターン1の人数が33名,火曜日からが26名...という結果になりました.
シフト計画(解)の表示
print("----- 結果 -----")
print("x1 : ", pulp.value(x1))
print("x2 : ", pulp.value(x2))
print("x3 : ", pulp.value(x3))
print("x4 : ", pulp.value(x4))
print("x5 : ", pulp.value(x5))
print("x6 : ", pulp.value(x6))
print("x7 : ", pulp.value(x7))
----- 結果 ----- x1 : 33.0 x2 : 26.0 x3 : 8.0 x4 : 26.0 x5 : 7.0 x6 : 27.0 x7 : 7.0
最小化された雇用すべき合計人数(目的関数の値)を表示します.これは上の結果の合計に一致するはずです.
目的関数の値
print("合計人数 = ", pulp.value(problem.objective))
合計人数 = 134.0
プログラム全体を示します.
プログラム全体
import pulp
import numpy as np
# 曜日ごとの必要人数の設定
d = np.array([100, 100, 100, 100, 100, 94, 75])
# シフトパターンの行列
A = np.array([
[1, 0, 0, 1, 1, 1, 1],
[1, 1, 0, 0, 1, 1, 1],
[1, 1, 1, 0, 0, 1, 1],
[1, 1, 1, 1, 0, 0, 1],
[1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1],
])
# 問題オブジェクトを作る
problem = pulp.LpProblem(name='シフト計画', sense=pulp.LpMinimize)
# 変数の定義
x1 = pulp.LpVariable('x1', lowBound=0, cat=pulp.LpInteger)
x2 = pulp.LpVariable('x2', lowBound=0, cat=pulp.LpInteger)
x3 = pulp.LpVariable('x3', lowBound=0, cat=pulp.LpInteger)
x4 = pulp.LpVariable('x4', lowBound=0, cat=pulp.LpInteger)
x5 = pulp.LpVariable('x5', lowBound=0, cat=pulp.LpInteger)
x6 = pulp.LpVariable('x6', lowBound=0, cat=pulp.LpInteger)
x7 = pulp.LpVariable('x7', lowBound=0, cat=pulp.LpInteger)
# 目的関数を設定する
problem += x1 + x2 + x3 + x4 + x5 + x6 + x7
for i in range(len(d)):
problem += A[i,0]*x1 + A[i,1]*x2 + A[i,2]*x3 + A[i,3]*x4 + A[i,4]*x5 + A[i,5]*x6 + A[i,6]*x7 >= d[i]
# 問題を表示する
print(problem)
# 問題を解く
#problem.solve()
problem.solve(pulp.PULP_CBC_CMD(msg=0)) # ログを出力しない
# 結果を表示する
print("Status : ", pulp.LpStatus[problem.status])
print("----- 結果 -----")
print("x1 : ", pulp.value(x1))
print("x2 : ", pulp.value(x2))
print("x3 : ", pulp.value(x3))
print("x4 : ", pulp.value(x4))
print("x5 : ", pulp.value(x5))
print("x6 : ", pulp.value(x6))
print("x7 : ", pulp.value(x7))
print("合計人数 = ", pulp.value(problem.objective))
シフト計画: MINIMIZE 1*x1 + 1*x2 + 1*x3 + 1*x4 + 1*x5 + 1*x6 + 1*x7 + 0 SUBJECT TO _C1: x1 + x4 + x5 + x6 + x7 >= 100 _C2: x1 + x2 + x5 + x6 + x7 >= 100 _C3: x1 + x2 + x3 + x6 + x7 >= 100 _C4: x1 + x2 + x3 + x4 + x7 >= 100 _C5: x1 + x2 + x3 + x4 + x5 >= 100 _C6: x2 + x3 + x4 + x5 + x6 >= 94 _C7: x3 + x4 + x5 + x6 + x7 >= 75 VARIABLES 0 <= x1 Integer 0 <= x2 Integer 0 <= x3 Integer 0 <= x4 Integer 0 <= x5 Integer 0 <= x6 Integer 0 <= x7 Integer Status : Optimal ----- 結果 ----- x1 : 33.0 x2 : 26.0 x3 : 8.0 x4 : 26.0 x5 : 7.0 x6 : 27.0 x7 : 7.0 合計人数 = 134.0
上のコードは変数や制約式を一つひとつ定義している関係で美しくなく,また拡張性に欠ける状態です.まずプログラムの構造を理解するという最初の段階としては上のような理解が必要ですが,より一般的にするためにはコードを改善すべきでしょう.
上のコードはリストや繰り返し,内包表記を使うと,より美しく拡張性のあるコードに書き直すことができます.最後には各曜日の出勤者数の合計を表示しています.水曜日以外の出勤者数は必要人数の下限と一致しており,水曜日は必要人数よりも1名多い作業員が出勤することが確認できました.(次の SciPy を使って解いた場合は日曜日が1名余分となるような別の結果が得られます.)
プログラム全体
import pulp
import numpy as np
# 曜日の一覧
yobi = np.array(["月", "火", "水", "木", "金", "土", "日"])
# 曜日ごとの必要人数の設定
d = np.array([100, 100, 100, 100, 100, 94, 75])
# シフトパターンの設定
A = np.array([
[1, 0, 0, 1, 1, 1, 1],
[1, 1, 0, 0, 1, 1, 1],
[1, 1, 1, 0, 0, 1, 1],
[1, 1, 1, 1, 0, 0, 1],
[1, 1, 1, 1, 1, 0, 0],
[0, 1, 1, 1, 1, 1, 0],
[0, 0, 1, 1, 1, 1, 1],
])
# 問題オブジェクトを作る
problem = pulp.LpProblem(name='シフト計画', sense=pulp.LpMinimize)
# 変数の定義
x = [pulp.LpVariable(f"x{i}", lowBound=0, cat=pulp.LpInteger) for i in range(1, A.shape[1]+1)]
# 目的関数を設定する
problem += sum(x)
# 制約式
for i in range(len(d)):
problem += pulp.lpSum(A[i][j] * x[j] for j in range(len(x))) >= d[i]
# 問題を表示する
print(problem)
# 問題を解く
# problem.solve()
problem.solve(pulp.PULP_CBC_CMD(msg=0)) # ログを出力しない
# 結果を表示する
print("Status : ", pulp.LpStatus[problem.status])
print("----- シフト計画の最適解 -----")
for xj in x:
print(f"{xj.name}: {xj.varValue}")
print("合計人数 = ", pulp.value(problem.objective))
print("----- 曜日ごとの出勤者数と必要人数 -----")
for i in range(len(d)):
n = pulp.lpSum(A[i][j] * x[j] for j in range(len(x))).value()
print(f"{yobi[i]} : {n} : {d[i]}")
シフト計画:
MINIMIZE
1*x1 + 1*x2 + 1*x3 + 1*x4 + 1*x5 + 1*x6 + 1*x7 + 0
SUBJECT TO
_C1: x1 + x4 + x5 + x6 + x7 >= 100
_C2: x1 + x2 + x5 + x6 + x7 >= 100
_C3: x1 + x2 + x3 + x6 + x7 >= 100
_C4: x1 + x2 + x3 + x4 + x7 >= 100
_C5: x1 + x2 + x3 + x4 + x5 >= 100
_C6: x2 + x3 + x4 + x5 + x6 >= 94
_C7: x3 + x4 + x5 + x6 + x7 >= 75
VARIABLES
0 <= x1 Integer
0 <= x2 Integer
0 <= x3 Integer
0 <= x4 Integer
0 <= x5 Integer
0 <= x6 Integer
0 <= x7 Integer
Status : Optimal
----- シフト計画の最適解 -----
x1: 33.0
x2: 26.0
x3: 8.0
x4: 26.0
x5: 7.0
x6: 27.0
x7: 7.0
合計人数 = 134.0
----- 曜日ごとの出勤者数と必要人数 -----
月 : 100.0 : 100
火 : 100.0 : 100
水 : 101.0 : 100
木 : 100.0 : 100
金 : 100.0 : 100
土 : 94.0 : 94
日 : 75.0 : 75
シフト計画問題は線形最適化問題とは異なり,決定変数(の一部/全体)が整数値しか取れないため,混合整数線形最適化問題になります.SciPyではこのような問題を scipy.optimize.milp
関数を使って解くことができます.この scipy.optimize.milp
関数で解くことができる問題の形式は
という形式です.この形式に従ってシフト計画問題を記述すると,目的関数では \({\bf{c}}^{\rm{T}} = (1, 1, 1, 1, 1, 1, 1)\) とすると良いでしょう.一方で,制約式の下限 \({\bf{d_l}}\) が出勤者数の下限となり,上限 \({\bf{d_u}}\) は存在しない,あるいは上限が無限大であると考えればよいでしょう.
まずは必要なモジュールをインポートします.
モジュールのインポート
import numpy as np
from scipy.optimize import milp
from scipy.optimize import LinearConstraint
結果表示で曜日がわかるように曜日のNumpy配列を準備します.
曜日のリスト
yobi = np.array(["月", "火", "水", "木", "金", "土", "日"])
目的関数の係数ベクトル \({\bf{c}}^{\rm{T}} = (1, 1, 1, 1, 1, 1, 1)\) を設定します.
目的関数の係数
c = np.array([1,1,1,1,1,1,1])
シフトパターンの行列 \(A\) を定義します.
シフトパターン
A = np.array([
[1,0,0,1,1,1,1],
[1,1,0,0,1,1,1],
[1,1,1,0,0,1,1],
[1,1,1,1,0,0,1],
[1,1,1,1,1,0,0],
[0,1,1,1,1,1,0],
[0,0,1,1,1,1,1]
])
制約条件の下限をベクトル形式で定義します.出勤者数の下限を設定すると良いことに注意してください.
制約条件の下限
d_l = np.array([100,100,100,100,100,94,75])
次に制約条件の上限をベクトル形式で定義します.ここでは,無限大のベクトルを定義しますが,一旦ゼロベクトルを作成してから,すべての要素を無限大に書き換えていることに注意してください.
制約条件の上限
d_u = np.zeros(len(d_l)) # 一旦 0 で初期化
d_u[:] = np.inf # すべての要素を np.inf に書き換え
制約条件を定義します.線形制約であるので,LinearConstraint
を利用し,下限 (Lower Bound) と上限 (Upper Bound) をそれぞれ指定します.
制約条件
constraints = LinearConstraint(A, lb=d_l, ub=d_u)
次に決定変数の整数制約を設定します.今回の問題では決定変数(各曜日の雇用者数)は整数であるので,すべて1を設定します.(0を指定すると実数の変数になります.詳細はこちらで確認してください.)
変数条件
integrality = np.array([1, 1, 1, 1, 1, 1, 1])
milp
関数を使ってシフト計画問題を解きます.Pulpで解いた結果と x の値が僅かに異なりますが,総数で134人を雇用するという結果は同じです.
問題を解く
result = milp(c, constraints=constraints, integrality=integrality)
print(result)
message: Optimization terminated successfully. (HiGHS Status 7: Optimal) success: True status: 0 fun: 134.0 x: [ 3.200e+01 2.600e+01 8.000e+00 2.600e+01 8.000e+00 2.600e+01 8.000e+00] mip_node_count: 1 mip_dual_bound: 134.0 mip_gap: 0.0
結果からシフト計画の最適解を取り出して表示します.なお,決定変数が x0
から x6
となっていることに注意してください.
最適解
print("----- シフト計画の最適解 -----")
for i, x in enumerate(result.x):
print(f"x{i}: {x}")
print("合計人数 = ", result.fun)
----- シフト計画の最適解 ----- x0: 32.0 x1: 26.0 x2: 8.0 x3: 26.0 x4: 8.0 x5: 26.0 x6: 8.0 合計人数 = 134.0
曜日ごとの出勤者数と必要人数を表示します.日曜日だけ1名余分に出勤し,その他の曜日については最低人数と一致していることがわかります.(Pulpで解いた場合は水曜日だけが1名余分でした.どちらも最適解であることは理解できると思います.)
曜日ごとの結果
print("----- 曜日ごとの出勤者数と必要人数 -----")
for i in range(len(d_u)):
n = np.sum([A[i][j] * result.x[j] for j in range(len(result.x))])
print(f"{yobi[i]} : {n} : {d_l[i]}")
----- 曜日ごとの出勤者数と必要人数 -----
月 : 100.0 : 100
火 : 100.0 : 100
水 : 100.0 : 100
木 : 100.0 : 100
金 : 100.0 : 100
土 : 94.0 : 94
日 : 76.0 : 75
なお,変数制約で0を指定した場合は実数の変数制約になります.この状態で問題を解くと次のような結果になることにも注意してください.つまり,シフトパターン1の雇用者数が 32.4名,シフトパターン2が 26.4名,...,となり,合計の雇用者数が 133.8名です.今回のシフト計画問題では実数値は認められないので,変数の整数制約を付与する必要があり,合計の134名が最適解であることも確認してください.
実数の変数制約
integrality = np.array([0, 0, 0, 0, 0, 0, 0])
# integrality = np.array([1, 1, 1, 1, 1, 1, 1])
result = milp(c, constraints=constraints, integrality=integrality)
print(result)
message: Optimization terminated successfully. (HiGHS Status 7: Optimal) success: True status: 0 fun: 133.8 x: [ 3.240e+01 2.640e+01 7.400e+00 2.640e+01 7.400e+00 2.640e+01 7.400e+00] mip_node_count: None mip_dual_bound: None mip_gap: None
ここで,SciPyを利用したシフト計画問題のプログラム全体を示します.
プログラム全体
import numpy as np
from scipy.optimize import milp
from scipy.optimize import LinearConstraint
# 曜日の一覧
yobi = np.array(["月", "火", "水", "木", "金", "土", "日"])
# 目的関数の係数
c = np.array([1,1,1,1,1,1,1])
# シフトパターン
A = np.array([
[1,0,0,1,1,1,1],
[1,1,0,0,1,1,1],
[1,1,1,0,0,1,1],
[1,1,1,1,0,0,1],
[1,1,1,1,1,0,0],
[0,1,1,1,1,1,0],
[0,0,1,1,1,1,1]
])
# 下限
d_l = np.array([100,100,100,100,100,94,75])
# 上限
d_u = np.zeros(len(d_l)) # 一旦 0 で初期化
d_u[:] = np.inf # すべての要素を np.inf に書き換え
# 線形制約を定義
constraints = LinearConstraint(A, lb=d_l, ub=d_u)
# 変数の整数制約(0は実数,1は整数)
# integrality = np.array([0, 0, 0, 0, 0, 0, 0]) # 実数制約の場合
integrality = np.array([1, 1, 1, 1, 1, 1, 1]) # 整数制約(今回はこれ)
# integrality = np.array([0, 0, 1, 1, 1, 0, 1]) # 実数と整数の混合問題の場合
# 問題を解く (Mixed-integer linear programming)
result = milp(c, constraints=constraints, integrality=integrality)
# 結果を表示する
print("\n----- 結果の表示 -----")
print(result)
# print("message : ", result.message)
print("----- シフト計画の最適解 -----")
for i, x in enumerate(result.x):
print(f"x{i}: {x}")
print("合計人数 = ", result.fun)
print("----- 曜日ごとの出勤者数と必要人数 -----")
for i in range(len(d_u)):
n = np.sum([A[i][j] * result.x[j] for j in range(len(result.x))])
print(f"{yobi[i]} : {n} : {d_l[i]}")
----- 結果の表示 ----- message: Optimization terminated successfully. (HiGHS Status 7: Optimal) success: True status: 0 fun: 134.0 x: [ 3.200e+01 2.600e+01 8.000e+00 2.600e+01 8.000e+00 2.600e+01 8.000e+00] mip_node_count: 1 mip_dual_bound: 134.0 mip_gap: 0.0 ----- シフト計画の最適解 ----- x0: 32.0 x1: 26.0 x2: 8.0 x3: 26.0 x4: 8.0 x5: 26.0 x6: 8.0 合計人数 = 134.0 ----- 曜日ごとの出勤者数と必要人数 ----- 月 : 100.0 : 100 火 : 100.0 : 100 水 : 100.0 : 100 木 : 100.0 : 100 金 : 100.0 : 100 土 : 94.0 : 94 日 : 76.0 : 75