定式化とソルバー

準備

from gurobipy import Model, quicksum, GRB
# from mypulp import Model, quicksum, GRB
from collections import OrderedDict, defaultdict
import networkx as nx

シフト最適化問題

様々な最適化のベンチマーク問題例を収集しているサイトとしては、昔からあるOR-Library (http://people.brunel.ac.uk/~mastjjb/jeb/info.html) が有名だ。 1990年に開始したこのライブラリには、論文で実験された問題のデータが配布されている。昔は、専用のアルゴリズムを作成してやっと解けていた問題が、 今では数理最適化ソルバーを使って比較的簡単に解けるようになっている。 簡単な練習問題ではなく、もっと難しい問題に挑戦したい人には、このライブラリはうってつけだ。

ここでは、以下のシフト最適化問題を解いてみる.

複数のジョブを作業員に割り当てたい。ジョブには開始時刻と終了時刻が与えられており、ジョブの処理には1人の作業員が必要で、 同時に2つのジョブの処理はできないものとする。作業員には固定費用があり、使用する作業員の総固定費用を最小化せよ。

もとになった論文やデータは以下ののサイトを参照されたい。

http://people.brunel.ac.uk/~mastjjb/jeb/orlib/ptaskinfo.html

定式化

  • $J$ : ジョブの集合
  • $W$: 作業員の集合
  • $b_w$: 作業員 $w$ の固定費用(ベンチマーク問題例ではすべて $1$)
  • $s_j, f_j$: ジョブ $j$ の開始時刻と終了時刻
  • $J_w$: 作業員 $w$ が遂行できるジョブの集合
  • $W_j$: ジョブ $j$ を処理できる作業員の集合
  • $G=(J,E)$: ジョブを点とし、2つの区間 $[s_i,f_i]$ と $[s_j,f_j]$ に共通部分があったときに、枝 $(i,j)$ をはったグラフ(区間グラフ)
  • $C$ : 区間グラフ $G=(J,E)$ に対する極大クリークの集合
  • $J_c$: クリーク(完全部分グラフ) $c$ に含まれるジョブの集合
  • $x_{jw}$: ジョブ $j$ が作業員 $w$ に割り当てられるとき1, それ以外のとき 0
  • $y_w$: 作業員 $w$ が使われるとき1, それ以外のとき0
$$ \begin{array}{lll} minimize & \sum_{w \in W} b_w y_y & \\ s.t. & \sum_{w \in W_j} x_{jw} =1 & \forall j \in J \\ & \sum_{j \in (J_c \cap J_w)} x_{jw} \leq y_w & \forall w \in W, c \in C \\ \end{array} $$

最初の制約は、すべてのジョブが処理されなければならないことを表す。 2番目の制約は、クリークに含まれるジョブは、高々1人の作業員にしか割り当てられないことと割り当てられた作業員は使われなければならないことを同時に表す。

まず,データを読み込んで準備をする.

fname = "../data/ptask/data_10_51_111_66.dat"
f = open(fname, "r")
lines = f.readlines()
f.close()
n_jobs = int(lines[4].split()[2])
print("number of jobs=", n_jobs)
start, finish = OrderedDict(), OrderedDict()
for j in range(n_jobs):
    start[j], finish[j] = map(int, lines[5 + j].split())
number of jobs= 111
n_workers = int(lines[5 + n_jobs].split()[2])
print("number of workers=", n_workers)
Jobs = defaultdict(set)
for w in range(n_workers):
    line = lines[6 + n_jobs + w].split()
    for j in line[1:]:
        Jobs[w].add(int(j))
number of workers= 51

2つのジョブの開始時刻と終了時刻に重なりがあるかどうか判定する関数を作っておく.

def intersect(j, k):
    if start[j] > start[k]:
        j, k = k, j
    if finish[j] > start[k]:
        return True
    else:
        return False

ジョブを点とした交差グラフを生成する.このグラフのクリーク(完全部分グラフ)が,同時に処理できないジョブの集合になる.

G = nx.Graph()
for j in range(n_jobs):
    for k in range(j + 1, n_jobs):
        if intersect(j, k):
            G.add_edge(j, k)

networkXを用いて極大クリークを列挙しておく.

Cliques = []
for c in nx.find_cliques(G):
    Cliques.append(set(c))
print("Number of cliques=", len(Cliques))
Number of cliques= 32
model = Model()
x, y = {}, {}
for w in range(n_workers):
    y[w] = model.addVar(vtype="B", name=f"y[{w}]")
    for j in Jobs[w]:
        x[j, w] = model.addVar(vtype="B", name=f"x[{j},{w}]")
model.update()
for j in range(n_jobs):
    model.addConstr(
        quicksum(x[j, w] for w in range(n_workers) if (j, w) in x) == 1,
        name=f"job_assign[{j}]",
    )
for w in range(n_workers):
    for j, c in enumerate(Cliques):
        model.addConstr(
            quicksum(x[j, w] for j in Jobs[w].intersection(c)) <= y[w],
            name=f"connection[{j},{w}]",
        )
model.setObjective(quicksum(y[w] for w in range(n_workers)), GRB.MINIMIZE)
model.optimize()
Academic license - for non-commercial use only - expires 2023-06-24
Using license file /Users/mikiokubo/gurobi.lic
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 1743 rows, 3943 columns and 48251 nonzeros
Model fingerprint: 0x6315480a
Variable types: 0 continuous, 3943 integer (3943 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 51.0000000
Presolve removed 449 rows and 0 columns
Presolve time: 0.08s
Presolved: 1294 rows, 3943 columns, 36453 nonzeros
Variable types: 0 continuous, 3943 integer (3943 binary)

Root relaxation: objective 4.000000e+01, 2652 iterations, 0.17 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0      40.0000000   40.00000  0.00%     -    0s

Explored 0 nodes (3964 simplex iterations) in 0.36 seconds
Thread count was 16 (of 16 available processors)

Solution count 2: 40 51 

Optimal solution found (tolerance 1.00e-04)
Best objective 4.000000000000e+01, best bound 4.000000000000e+01, gap 0.0000%
print("Optimal value=", model.ObjVal)
Optimal value= 40.0

看護婦スケジューリング問題

典型的な看護婦スケジューリング問題(nurse scheduling problem)では,以下の制約が付加される.

  • 毎日の各勤務(昼、夕、夜)の必要人数
  • 各看護婦に対して30日間の各勤務日数の上下限
  • 指定休日、指定会議日
  • 連続7日間に、最低1日の休日、最低1日の昼
  • 禁止パターン:3連続夜、4連続夕、5連続昼、夜勤明けの休日以外、夕の直後の昼あるいは会議、休日--勤務--休日
  • 夜は2回連続で行う
  • 2つのチームの人数をできるだけ均等化

このような付加条件が課された看護婦スケジューリング問題は,制約最適化ソルバーSCOPで簡単に定式化できる.

看護婦 $i$ の 日 $d$におけるシフトを表す変数 $X_{id}$ (領域は昼,夕,夜,休暇の集合)を用いる.

対応する値変数 $x_{ids}$ は $0$-$1$ 変数であり, 看護婦 $i$ の $d$ 日のシフトが $s$ のとき $1$ になる.

制約を式で表現する.

  • 毎日の各勤務(昼、夕、夜)の必要人数

日 $d$ のシフト $s$ の必要人数を $L_{ds}$ とする. $$ \sum_{i} x_{ids} \geq L_{ds} \ \ \ \forall d,s $$

  • 各看護婦に対して30日間の各勤務日数の上下限

休暇以外のシフト集合を $W$,下限と上限を $LB,UB$ とする. $$ LB \leq \sum_{d, s \in W} x_{ids} \leq UB \ \ \ \forall i $$

  • 指定休日、指定会議日

看護婦 $i$ が休日を希望する日の集合を $R_i$ とする. $$ \sum_{d \in R_i, s \in W} x_{ids} \leq 0 \ \ \ \forall i $$

  • 連続7日間に、最低1日の休日、最低1日の昼

$d$ 日から始まる連続7日間の集合を $C7_d$ とする.休日(値は $r$) の制約だけ記述する. $$ \sum_{t \in C7_d} x_{itr} \geq 1 \ \ \ \forall i, d $$

  • 禁止パターン:3連続夜、4連続夕、5連続昼、夜勤明けの休日以外、夕の直後の昼あるいは会議、休日-勤務-休日

例として3連続夜勤の場合の制約を考える.夜勤の値を $n$, $d$ 日から始まる連続3日間の集合を $C3_d$ とする.

$$ \sum_{t \in C3_d} x_{itn} \leq 2 \ \ \ \forall i, d $$
  • 夜は2回連続で行う

夜勤以外のシフトの集合を $N$ とする. $$ \sum_{s \in N} x_{ids} + x_{i,d+1,n} + \sum_{s \in N} x_{i,d+2,s} \leq 2 \ \ \ \forall i, d $$

  • 2つのチームの人数をできるだけ均等化

1つのチームに含まれる看護婦の集合を $T$ とし,以下の制約を考慮制約とする. $$ \sum_{i \in T} x_{ids} = \sum_{i \not\in T} x_{ids} \ \ \ \forall d, s $$

  • シフトの開始期と終了期の情報がほしい場合

以下の変数を追加する.

  • $Z_{id}$: 看護婦 $i$ が日 $d$ に開始するシフト(連続してシフトを行う場合のダミーのシフトを領域に追加しておく). 対応する値変数 $z_{ids}$ は,看護婦 $i$ がシフト $s$ を日 $d$ に開始するとき $1$

  • $W_{id}$: 看護婦 $i$ が日 $d$ に終了するシフト(連続してシフトを行う場合のダミーのシフトを領域に追加しておく). 対応する値変数 $w_{ids}$は,看護婦 $i$ がシフト $s$ を日 $d$ に終了するとき $1$

$$ \begin{array}{l l } z_{ids} \leq x_{ids} & \forall i,d,s \\ z_{ids} - w_{i,d-1,s} = x_{ids} -x_{i,d-1,s} & \forall i,d,s \end{array} $$

看護婦スケジューリングの国際コンペティション(First International Nurse Rostering Competition 2010)では, SCOPは3つの異なるトラックですべてメダル(2位,3位,4位)を獲得している.

業務割り当てを考慮したシフトスケジューリング問題

実際のシフト最適化においては,複数の日(たとえば1ヶ月)にスタッフを割り当てる必要がある.このとき,各日の各時間帯においてどのような業務をどのスタッフに割り当てるかを決めたい. もちろん,日や時間帯によって,各業務の遂行に必要な人数は異なる.このような諸条件を考慮した最適化は,問題依存であり,一般には難しい.

実務から発生した諸条件を組み込んだシフトシフト最適化システムとして OptShift (https://www.logopt.com/optshift/) が開発されている. これは,制約最適化ソルバー SCOP を用いたものであり,大規模な実際問題に対する実用的な解を,現実的な計算時間で算出する.