ある工場では平日に100名以上,土曜日に90名以上,日曜日に80名以上の作業員を必要としている.各作業員は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 というモジュールを使えば,このようなシフト計画問題を Python で簡単かつ高速に解くことができます.
PuLP は Anaconda をインストールしただけではインストールされていないので,まずは PuLP をインストールします.PuLP のインストールは「Anaconda Prompt」または「Anaconda Powershell Prompt」で pip install pulp
を実行すれば良いでしょう.なお,学内の情報処理実習室では pip install pulp --user
のように --user
オプションを指定する必要があるかもしれません.また,インストールされているモジュールの一覧は pip list
で確認できます.Jupyter Notebook のマジックコマンドを利用しても構いません.
まずは一つずつ変数や制約式を定義して PuLP でシフト計画問題を解いてみよう.
PuLP モジュールをインポートします.もしもエラーが表示されたら,準備を参考に,PuLP をインストールしよう.また,行列を利用するので,NumPyもインポートします.
モジュールのインポートimport pulp
import numpy as np
曜日ごとに必要な人数を格納するためのベクトルを定義します.
曜日ごとの必要人数の設定d = np.array([100, 100, 100, 100, 100, 90, 80])
シフトパターンの行列を定義します.
シフトパターンの行列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 >= 90 _C7: x3 + x4 + x5 + x6 + x7 >= 80 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.PULP_CBC_CMD(msg=0)) # ログを出力しない
1
問題を解けたかどうか念のために確認します.Optimal
と表示されたので正しく解けているようです.
結果の確認print("Status : ", pulp.LpStatus[problem.status])
Status : Optimal
最後に結果を表示します.月曜日から出勤するシフトパターン1の人数が32名,火曜日からが22名...という結果になりました.
シフト計画(解)の表示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 : 32.0 x2 : 22.0 x3 : 12.0 x4 : 22.0 x5 : 12.0 x6 : 22.0 x7 : 12.0
最小化された雇用すべき合計人数(目的関数の値)を表示します.これは上の結果の合計に一致するはずです.
目的関数の値print("合計人数 = ", pulp.value(problem.objective))
合計人数 = 134.0
プログラム全体を示します.
プログラム全体import pulp
import numpy as np
# 曜日ごとの必要人数の設定
d = np.array([100, 100, 100, 100, 100, 90, 80])
# シフトパターンの行列
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.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 >= 90 _C7: x3 + x4 + x5 + x6 + x7 >= 80 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 : 32.0 x2 : 22.0 x3 : 12.0 x4 : 22.0 x5 : 12.0 x6 : 22.0 x7 : 12.0 合計人数 = 134.0
上のコードは変数や制約式を一つひとつ定義している関係で美しくなく,また拡張性に欠ける状態です.まずプログラムの構造を理解するという最初の段階としては上のような理解が必要ですが,より一般的にするためにはコードを改善すべきでしょう.
上のコードはリストや繰り返し,内包表記を使うと,より美しく拡張性のあるコードに書き直すことができます.最後には各曜日の出勤者数の合計を表示しています.これらの値は必要人数の下限と一致していることも確認できました.
プログラム全体import pulp
import numpy as np
# 曜日の一覧
yobi = np.array(["月", "火", "水", "木", "金", "土", "日"])
# 曜日ごとの必要人数の設定
d = np.array([100, 100, 100, 100, 100, 90, 80])
# シフトパターンの設定
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.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.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 >= 90 _C7: x3 + x4 + x5 + x6 + x7 >= 80 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: 32.0 x2: 22.0 x3: 12.0 x4: 22.0 x5: 12.0 x6: 22.0 x7: 12.0 合計人数 = 134.0 ----- 曜日ごとの出勤者数と必要人数 ----- 月 : 100.0 : 100 火 : 100.0 : 100 水 : 100.0 : 100 木 : 100.0 : 100 金 : 100.0 : 100 土 : 90.0 : 90 日 : 80.0 : 80