組合せ最適化への挑戦録（勤務スケジューリング）6.びっくりするほど回り道

はじめに

なんか色々とやろうとして失敗した記録なので読み飛ばしても問題ないやつです

これまでの踏跡

目次：今回の取組課題

目的関数の修正プラン
参考サイトの理解
目的関数の修正
補足：プログラム全容

目的関数の修正プラン

現状、勤務スケジューリングで必要な目的関数は下記の4つである。
1. 合計勤務回数はなるべく少なくする方が良い

1. 勤務者間の勤務回数はなるべく等しいほうが良い

2. 有給の希望がなるべく通る方が良い

3. 夜勤忌避の希望がなるべく通る方が良い

1.3.4.については、おおよその見当がついており、

for d,a in list_da:
  obj1_temp = pulp.lpSum(x[p,m,d,a] for p,m in list_pm)
  obj1=obj1+obj1_temp

上記のように条件を満たす変数の総和を引っ張ってきて、最終的に、

ShiftScheduling += obj1+obj2+obj3+obj4,"Target"

と、目的変数がobj1~4を最小にすることですよ、と入れてやれば成立するはずである。
後は優先順位として、1.3.4に定数倍して重み付けをする。

問題は2. 勤務者間の勤務回数はなるべく等しいほうが良い、である。
普通に考えると、分散か標準偏差あたりを最小化すれば良いのだが……

参考サイト:
1 https://qiita.com/SaitoTsutomu/items/f4478dfbc3c1cf6425e3
2 https://qiita.com/matsulib/items/898873b73d584c7dcb8b

上記がかなり難解で、すぐに理解できなかったので、これの勉強にかなり時間を要した。
せっかくなので纏めておく。
なお、参考サイト1にある平均偏差最小化は未だに良く分かっていないので、ここでは参考サイト2の範囲最小化についてまとめる。

参考サイトの理解

参考サイト2から引用

m = LpProblem() # 数理モデル
x = np.array(addbinvars(メンバ数,チーム数)) # 割当
y = np.array(addvars(チーム数,2)) # チーム内の最小と最大
z = addvars(2) # チーム間の最小と最大
m += lpSum(y[:,1]-y[:,0]) + 1.5*(z[1]-z[0]) # 目的関数
for i in range(メンバ数):
    m += lpSum(x[i]) == 1 # どこかのチームに所属

---こっからよく分からんゾーン---
for j in range(チーム数):
    m += lpDot(得点.sum(1),x[:,j]) >= z[0]
    m += lpDot(得点.sum(1),x[:,j]) <= z[1]
    for k in range(得点.shape[1]):
        m += lpDot(得点.iloc[:,k],x[:,j]) >= y[j,0]
        m += lpDot(得点.iloc[:,k],x[:,j]) <= y[j,1]

まず前半部は多少なり理解できる。
xは6行3列の変数行列、yは3行2列の変数行列、zが2個の変数配列、というのは分かる。

目的関数が「yの1列目-0列目」と「zの1列目-0列目」の最小化というのも分かる。
ただ、これ以降の数式を見て「yとzの1列目にはチームの最大スコア、0列目には最小スコアが入る」というのが分からなかった。

可能な限り分解して捉えていく。
まずは各変数とdataframeの中身を確認すると、

このようになっている。
実際の変数xyzには行ラベル、列ラベルはないが、ここでは定義的な意味合いを確認するために明記した。

では、実際にfor文を追ってみよう。

for j in range(チーム数):
    m += lpDot(得点.sum(1),x[:,j]) >= z[0]
    m += lpDot(得点.sum(1),x[:,j]) <= z[1]

jはrange()関数なので、[0,1,2]とカウントアップされるだけのシンプルなものだ。
で、問題の2行目。

lpDotは行列同士の内積を表す。
得点.sum(1)は得点を行ごとに加算した結果である。
x[ : , j ]は変数xのj行目。

j=1を代入した場合、

このようになる。
この時、Sum(Co...An)は各メンバーの得点総和を表し、xは「メンバーがそのチームに所属しているかどうか？」を表すため、

lpDot(得点.sum(1),x[:,j])　は各チームの持つ総得点を表す。
この時、チーム内のメンバーの組み合わせは未確定である

ということになる。

for j in range(チーム数):
    m += lpDot(得点.sum(1),x[:,j]) >= z[0]
    m += lpDot(得点.sum(1),x[:,j]) <= z[1]

そしていよいよ本題の「>=z[0]」と「<=z[1]」となる。

普通に読めば、この式は
各チーム（ABC)の持つ総得点（P_A,P_B,P_C）は、全てz[0]以上、z[1]以下。
　という意味に他ならない。

そのため、この式からz[0]はチーム得点の最小値、z[1]はチーム得点の最大値が
　導出されて然るべきなのだが・・・いまいち腑に落ちない。

ここの概念がいまいち掴めておらず、
「未確定の点数(P_A~C)を用いて、未確定の変数(Z[0,1]）を定義する」
という式に対して、何でこれで決めていいの？感が凄い。

変数の形式（非負の連続変数など）を先に定義し、目的関数に放り込んだ後、
制約関数の式の中で、未確定の変数の総和から変数を定める・・・やっている事は確かにこの通りなのだが、
自分で立式せよ、と言われたときには全くピンとこない。

……非常に情けない話だが、3日ぐらい置いてみれば、ふと理解できたりすることもあるため、
　ひとまず丸パクリで作ろうと思う。

目的関数の修正

さて、ひとまず丸パクリしてみた結果が下記なのだが・・・解けない。
google colabの計算が何時まで経っても終わらない。
はてさて、どうしたものやら。
とりあえず今回は七転八倒した記録だけ残しておくことにする

プログラム全容

import openpyxl
import pandas as pd
import numpy as np
import pulp
from ortoolpy import addvars, addbinvars

#問題の定義
ShiftScheduling = pulp.LpProblem("ShiftScheduling", pulp.LpMinimize)

#変数宣言
#xは出欠を表し、13行(7日2シフトで7日夜だけ無し)、8列（勤務者8人）の0-1変数リスト
x = np.array(addbinvars(13,8))

#zは各勤務者の出勤回数の最大値と最小値を表す2個の変数リスト
z = addvars(2)

#目的関数の定義　全出席回数と、勤務者ごとの出勤回数の最大値-最小値、の合計を最小化する
ShiftScheduling += pulp.lpSum(x)+(z[1]-z[0]) ,"Target"

#制約条件1:全員必ず3回以上出勤
for m in range(8):
  ShiftScheduling += pulp.lpSum(x[:,m]) >= 3

#制約条件2 z[1]は勤務者の出勤回数の最大値、z[0]は出勤回数の最小値
for m in range(8):
  ShiftScheduling += pulp.lpSum(x[:,m]) >= z[0]
  ShiftScheduling += pulp.lpSum(x[:,m]) <= z[1]

#制約条件3 全シフトで2人以上出勤
for d in range(13):
  ShiftScheduling += pulp.lpSum(x[d,:]) >= 2

results = ShiftScheduling.solve()
print("optimality = {:}, target value = {:}".format(pulp.LpStatus[results], pulp.value(ShiftScheduling.objective)))

kekka=[]
for v in ShiftScheduling.variables():
#  print(f"{v} = {pulp.value(v)}")
  kekka = kekka+[int(pulp.value(v))]
print(kekka[-2:])
kekka=kekka[:-2]
kekka=pd.DataFrame(np.array(kekka).reshape(8,13).T)

print(kekka)

上記プログラムを動かすと、solve()のところで終わらなくなる。

考察

プログラムが動く条件を探してあちこち数字をいじっていくと、目的関数の重み付けを
大きく落としたら動いた。

ShiftScheduling += pulp.lpSum(x)/26+(z[1]-z[0]) ,"Target"

しかし・・・これで得られた結果がおかしい。

>>>出力結果
z=[4, 4]
    0  1  2  3  4  5  6  7
0   0  1  0  0  0  1  0  0
1   1  1  1  0  0  0  0  1
2   0  0  0  0  0  1  0  0
3   0  1  0  1  0  0  0  1
4   0  0  0  0  1  0  0  0
5   1  0  0  1  0  1  0  1
6   0  0  0  0  0  1  1  1
7   0  0  0  0  0  0  0  1
8   0  1  1  0  0  0  0  0
9   0  0  0  0  0  0  1  0
10  0  1  0  0  0  1  0  0
11  0  1  0  1  1  0  0  0
12  1  1  1  1  0  1  0  0

z:勤務者ごとの最小最大ともに4となっているが、出力表だと作業者1は明らかに7回働いていて、z[1]=4というのはおかしい。
袋小路に入ってしまった感じがある・・・7.ではもう一度練り直しから始める。進まねえなあ、おい