LoginSignup
1
2

More than 1 year has passed since last update.

@ookamikujira(たかひろ)

Pulp 組み合わせ最適化でシフト表作成　可能な限り2連休を作る

Last updated at Posted at 2022-05-27

注意
https://qiita.com/ki073/items/d45cd648e0625c488ee6
この記事よりもいい記事を見つけたのでこの記事を見ずに
どうぞ上記記事を参考にしてください。
GLPKを利用されてますが
この記事を見るより参考になりそうです。
自分もこの記事を参考に書きなおす予定です。

＜追記＞
参考にしてコードを書き直したら
コードもすっきしただけでなく
計算時間が短くなりました。

Python　Pulp での記述の仕方を下の方で紹介しています。
自分で考えたコードとの実行時間の比較もしております。
気になる方はご確認ください。

自分で考えたコード

可能な限り「１」を2連続で作りたい。

PythonとPulpを利用すればシフト表を自動で作成できますが(過去記事参照)
2連休で休みをできるだけ作りたいなんて時もありますよね。

シフト表を出勤を０と休みを１で表現すると
2連休は「１」「１」。
2連休を増やしたい。
3連休ではなく。

そんな人のため記事を書きました。

縦が日付、横が人として
つまり縦にに連続で１を作るイメージです。

最初に出力結果を見てもらえば何がしたいか分かりやすいかも。

コード

ちなみに少しでも大きな規模で計算させると計算時間が爆伸びします。
この記事で紹介する制約条件のみでは試していませんが
かなり複雑なシフト表を作成するコードに一緒に組み込んだ際には
30×23のマスで1時間2時間かけても最適解は見つからず。
（暫定最適解は出る。）

!pip install pulp

!pip install ortoolpy

import pandas as pd
from pulp import LpProblem , LpMinimize , LpVariable , lpSum , PULP_CBC_CMD , LpStatus
from ortoolpy import addvar, addvars, addbinvars

#モデルの定義/変数の定義
H = [i for i in range(12)]
W = [j for j in range(4)]
prob = LpProblem('TEST',LpMinimize)

X = pd.DataFrame(addbinvars(len(H), len(W)), index=H, columns=W)

#1の数を固定
for j in W:
  prob += lpSum(X.iloc[i,j] for i in H) == 6

#2連続「１」にインセンティブを与える
D = {(i,j) for i in H[:-1] for j in W}
ABS_one_two = LpVariable.dicts('ABS_one_two',D,cat='Binary') #0または1
Point_Duble = LpVariable.dicts('Point_Duble',D,lowBound=-2,cat='Integer') #-2または0

for i in H[:-1]:
  for j in W:
    prob += X.iloc[i,j] - X.iloc[i+1,j] <= ABS_one_two[i,j]
    prob += X.iloc[i,j] - X.iloc[i+1,j] >= -ABS_one_two[i,j]
    prob += ABS_one_two[i,j] - (X.iloc[i,j] + X.iloc[i+1,j] ) == Point_Duble[i,j]

#3連続「１」にペナルティを与える
T = {(i,j) for i in H[:-2] for j in W}
ABS_one_third = LpVariable.dicts('ABS_one_third',T,cat='Binary') #0または1
ABS_relay = LpVariable.dicts('ABS_relay',T,cat='Integer') #1または3
Triple_flg = LpVariable.dicts('Triple_flg',T,cat='Binary') #0または1

for i in H[:-2]:
  for j in W:
    prob += X.iloc[i,j] - X.iloc[i+2,j] <= ABS_one_third[i,j]
    prob += X.iloc[i,j] - X.iloc[i+2,j] >= -ABS_one_third[i,j]
    prob += 1 - (X.iloc[i,j] + 2*X.iloc[i+1,j] + X.iloc[i+2,j] - ABS_one_third[i,j]) <= ABS_relay[i,j]
    prob += 1 - (X.iloc[i,j] + 2*X.iloc[i+1,j] + X.iloc[i+2,j] - ABS_one_third[i,j]) >= -ABS_relay[i,j]
    prob += (ABS_relay[i,j] - 1)*4 == Triple_flg[i,j]*8



#目的関数(１が連続するたび-2点、3連続すると+8点)
prob.objective =  lpSum(Point_Duble[i,j]for i in H[:-1] for j in W)\
               +  8*lpSum(Triple_flg[i,j]for i in H[:-2] for j in W)

#計算
pulpCBC_setting = PULP_CBC_CMD(msg=1, threads=4, timeLimit=3600)
status = prob.solve(pulpCBC_setting)
print(LpStatus[status])
print(prob.objective.value())

#計算結果を成形・表示
new_shift = pd.DataFrame([[""for j in W]for i in H])
for i in H:
  for j in W:
      new_shift.iloc[i,j] = X.iloc[(i,j)].value()

new_shift_DF = pd.DataFrame(new_shift)
display(new_shift_DF)

出力

ちゃんと縦に2連続で「１」を量産できていますね。

Optimal
-24.0
    0	1	2	3
0	0.0	1.0	1.0	0.0
1	1.0	1.0	1.0	0.0
2	1.0	0.0	0.0	1.0
3	0.0	0.0	0.0	1.0
4	0.0	1.0	1.0	0.0
5	1.0	1.0	1.0	0.0
6	1.0	0.0	0.0	1.0
7	0.0	1.0	0.0	1.0
8	1.0	1.0	1.0	0.0
9	1.0	0.0	1.0	1.0
10	0.0	0.0	0.0	1.0
11	0.0	0.0	0.0	0.0

2連続「１」にインセンティブを与える

モデルがLpMinimizeとなっており最小化問題にしているため「１」が連続で現れたときに得点を下げるようにしてあげます。
結論を言ってしまうと、「１」が2連続で並ぶと-2点とします。

それぞれのマスには「０」または「１」が入ります。
2連続の「１」とは縦二つに並ぶマスに「１」が入っている状態です。

2連続のマスには
パターン１：「０」「０」
パターン２：「０」「１」
パターン３：「１」「０」
パターン４：「１」「１」
この4パターンです。

このパターン4の時だけマイナスポイントが付くようにしたいです。
これを実現するために
「マス1」と「マス2」の絶対値をとります。＝ ABS_one_two
次にABS_one_twoから「マス1」と「マス2」の和を引きます。
これを行うと
パターン１：「０」「０」→0
パターン２：「０」「１」→0
パターン３：「１」「０」→0
パターン４：「１」「１」→-2
パターン４が作られたとき-2点となります。

3連続「１」にペナルティを与える

上記のおかげで「１」が連続で並ぶようになりますが、これだけだと
2連続で止まらず6連続の「１」になります。

そのために3連続以上に対してペナルティを与えます。
もちろん3連続以上「１」を作成しない制約を設ければ早いのですが
現実のシフト表では3連続や4連続の希望休みが出たりします。
つまり希望があれば3連続や4連続も作りたいのです。

できるだけ3連続以上を作らないために
3連続の「１」を作った時にはペナルティを与えます。

結論を言ってしまうと3連続「１」につき+8点にしています。（絶対に8である必要はない。）

3連続「１」の中には2連続「１」が二つあるため
前述のインセンティブにより3連続「１」は-4点となっています。

そして2連続「１」が二つの場合も-4点です。3連続のほうがコスパがいい。
（利用する「１」が少なくてOKだから）
このため3連続「１」に対してはより大きなペナルティを与える必要があります。

ペナルティが8点であれば
インセンティブと合わせて
3連続「１」の場合、4点
2連続「１」の場合、-2点
それ以外の場合、0点
を作り出すことができます。

3つの連続するマスのパターンは以下8パターン。
パターン１：「０」「０」「０」
パターン２：「０」「０」「１」
パターン３：「０」「１」「０」
パターン４：「０」「１」「１」
パターン５：「１」「０」「０」
パターン６：「１」「０」「１」
パターン７：「１」「１」「０」
パターン８：「１」「１」「１」

このパターン８の時だけ+8点にしたいです。
ややこしいですがまず
「マス１」と「マス３」の差の絶対値をとります。＝　ABS_one_third
そして
1 - （「マス１」+　2*「マス２」+「マス３」- ABS_one_third ）　の絶対値をとります。＝　ABS_relay
最後に
(ABS_relay - 1)*4　をすると
パターン１：「０」「０」「０」→0
パターン２：「０」「０」「１」→0
パターン３：「０」「１」「０」→0
パターン４：「０」「１」「１」→0
パターン５：「１」「０」「０」→0
パターン６：「１」「０」「１」→0
パターン７：「１」「１」「０」→0
パターン８：「１」「１」「１」→8
となります。

多分もっとシンプルなやり方があると思うのですが
思いつきません。おしえて・・・

参考にしてコードの書き直し

出力される答えを同じにするため制約を一部追加しています。

####参考にして書き直したコード####
from ortoolpy import addvar, addvars, addbinvars
import pandas as pd
from pulp import LpProblem , LpMinimize , LpVariable , lpSum , PULP_CBC_CMD , LpStatus
import time

#モデルの定義/変数の定義
H = [i for i in range(10)]
W = [j for j in range(5)]
prob = LpProblem('TEST',LpMinimize)

X = pd.DataFrame(addbinvars(len(H), len(W)), index=H, columns=W)

#1の数を固定
for j in W:
  prob += lpSum(X.iloc[i,j] for i in H) == 6

#101を作れない
for i in H[:-2]:
  for j in W:
    prob +=X.iloc[i,j] + X.iloc[i+2,j] - X.iloc[i+1,j] <= 1

#2連続「１」を検出する
D = {(i,j) for i in H[:-1] for j in W}
Double_flg = LpVariable.dicts('Double_flg',D,cat='Binary') #0または1

for i in H[:-1]:
  for j in W:
    prob += X.iloc[i,j] + X.iloc[i+1,j] -1 <= Double_flg[i,j]
    prob += X.iloc[i,j] + X.iloc[i+1,j] >= 2*Double_flg[i,j]

#3連続「１」を検出する
T = {(i,j) for i in H[:-2] for j in W}
Triple_flg = LpVariable.dicts('Triple_flg',T,cat='Binary') #0または1

for i in H[:-2]:
  for j in W:
    prob += X.iloc[i,j] + X.iloc[i+1,j] + X.iloc[i+2,j] - 2 <= Triple_flg[i,j]
    prob += X.iloc[i,j] + X.iloc[i+1,j] + X.iloc[i+2,j] >= 3*Triple_flg[i,j]


#目的関数(１が連続するたび-2点、3連続すると+8点)
prob.objective =  -2*lpSum(Double_flg[i,j]for i in H[:-1] for j in W)\
               +  8*lpSum(Triple_flg[i,j]for i in H[:-2] for j in W)

#計測開始
start = time.time() 

#計算
pulpCBC_setting = PULP_CBC_CMD(msg=1, threads=4, timeLimit=3600)
status = prob.solve(pulpCBC_setting)
# 測定終了
elapsed = time.time() - start 
print(f"{elapsed * 1000:.1f} ms")  
print(LpStatus[status])
print(prob.objective.value())

#計算結果を成形・表示
new_shift = pd.DataFrame([[""for j in W]for i in H])
for i in H:
  for j in W:
      new_shift.iloc[i,j] = X.iloc[(i,j)].value()

new_shift_DF = pd.DataFrame(new_shift)
display(new_shift_DF)

出力

###参考にして書き直したコード結果
535.1 ms
Optimal
-30.0
　　0	1	2	3	4
0	1.0	1.0	1.0	1.0	1.0
1	1.0	1.0	1.0	1.0	1.0
2	0.0	0.0	0.0	0.0	0.0
3	0.0	0.0	0.0	0.0	0.0
4	1.0	1.0	1.0	1.0	1.0
5	1.0	1.0	1.0	1.0	1.0
6	0.0	0.0	0.0	0.0	0.0
7	0.0	0.0	0.0	0.0	0.0
8	1.0	1.0	1.0	1.0	1.0
9	1.0	1.0	1.0	1.0	1.0

###同じ条件で自分で書いたコード結果
1953.3 ms
Optimal
-30.0
　　0	1	2	3	4
0	1.0	1.0	1.0	1.0	1.0
1	1.0	1.0	1.0	1.0	1.0
2	0.0	0.0	0.0	0.0	0.0
3	0.0	0.0	0.0	0.0	0.0
4	1.0	1.0	1.0	1.0	1.0
5	1.0	1.0	1.0	1.0	1.0
6	0.0	0.0	0.0	0.0	0.0
7	0.0	0.0	0.0	0.0	0.0
8	1.0	1.0	1.0	1.0	1.0
9	1.0	1.0	1.0	1.0	1.0

出力された結果は同じですが
変数のシンプルさ、ゆえにコードも読みやすく
ロジックもわかりやすい。
そして計算時間が短くて済む。と完敗をしています。
(規模を大きくしたら計算時間はかなり伸びちゃうところは同じだけど)

連続「１」の抽出

X1 + X2 -1 <= flg
X1 + X2 >= 2*flg

X1 , X2 , flgはそれぞれ0か1の値をとるバイナリの変数です。
この二つの式が成り立つとき
X1 , X2 の論理積がflgとなります。
つまりX1 , X2の両方が1だった場合のみflgが1となります。

3連続も同じ要領です。

X1 + X2 + X3 - 2 <= flg
X1 + X2 + X3 >= 3*flg

###追記
二つの制約式を利用していましたがこれ一個の制約式でできると気が付いた。

X1 + X2 + X3 + x4 - 3 <= flg
X1 + X2 + X3 +x4 >= 4*flg

0110という並びがあればflgが立つ
このflgをの数を最大化すればよい。

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2