こんにちは。エンジニアの柴です。
今回は数理最適化問題を解いてみたいと思います。

数理最適化とは

与えられた制約条件の下で、現実の問題を数式で表し、目的関数の値を最大化（最小化）する解を1つ求めることです。

この数理最適化を用いると以下の有名な問題を解くことができます。

  • 巡回セールスマン問題
    セールスマンがある都市から出発し、すべての都市を一度ずつ訪れ、最初の都市に戻るルートのうち総移動距離が最小となるルートを求める問題。

  • ナップサック問題
    N個の品物があり、各品物に重さと価値が決められている。これらの品物から重さの総和がWを超えないようにいくつか選んだときの、価値の総和が最大となる選び方を求める問題。

また、身近な問題の例だと以下のような問題が挙げられます。

  • 配送計画問題
    最短の時間で宅配便を全員に配達するルートを求める問題

  • 職場での勤務シフト表の作成
    すべての従業員の勤務時間などの希望をみたしつつ勤務シフト表を作成する問題

数理最適化問題は連続型と離散型に大別され、上に挙げた有名な問題は離散最適化問題の整数計画問題に分類されます。
また、目的関数と制約条件がともに線形で表されるものを線形計画問題といいます。線形計画問題を簡単に言うと、連立方程式・不等式の変数が非常に多いときの解を求めるイメージです（100元の連立不等式・方程式を解くイメージです）。
整数計画問題と線形計画問題の両方の特徴を持つものを線形整数計画問題といいます。巡回セール万問題とナップザック問題もこの線形整数計画問題に含まれます。

線形計画問題の解を求める方法には厳密解法と近似解法に分けられます。厳密解法は唯一の最良の解を求める方法で、どうしても計算時間が長くなる傾向があります。これに対して近似解法は最良の解を求める保証はありませんが、計算時間が比較的短くなるという傾向があります。

今回は線形整数計画問題の勤務シフト表の作成を厳密解法で行いたいと思います。

条件の設定

今回は、シフトを作成する条件は以下のようにしたいと思います。
（職場に応じて条件を変更することが可能です）

  • 1日のシフトは早番、中番、遅番の3シフト
  • 早番は2人、中番は3人、遅番は2人が必要（過不足なくシフトに入れる）
  • 深夜勤務から早朝勤務にしないために、遅番→早番の順の勤務にならないようにする
  • 従業員のシフトと勤務回数の希望は以下の通り
    Aさん：早番、中番、遅番のどこでも入れる（週5～6回）
    Bさん：早番、中番、遅番のどこでも入れる（週4～5回）
    Cさん：早番、中番、遅番のどこでも入れる（週3～5回）
    Dさん：早番だけ入れる（週4～5回）
    Eさん：早番だけ入れる（週3～4回）
    Fさん：早番だけ入れる（週2～3回）
    Gさん：中番、遅番だけ入れる（週2～3回）
    Hさん：中番、遅番だけ入れる（週2～3回）
    Iさん：中番、遅番だけ入れる（週4～5回）
    Jさん：中番、遅番だけ入れる（週4～5回）
    Kさん：遅番だけ入れる（週4～5回）
    Lさん：遅番だけ入れる（週5～6回）

問題の定式化

線形計画問題では最適化（最大化または最小化）するための目的関数を定めます。
今回のシフト作成ではシフトの設定人数にぴったり人数を配置することになりますので、目的関数を0（最大値と最小値を0）に設定します。

定数の定義

$$ \begin{align} & \text{Shifts}_s \text{：各人の勤務可能なシフト（例：早番だけ入れるなど）} \\\\ & \text{People}_s \text{：シフト}s\text{に入れる人} \\\\ & \text{R}_s \text{：シフト}𝑠\text{に必要な人数（例：早番は2人など）} \end{align} $$

変数定義

制約条件で利用される変数を以下のように定義します。

$$ \begin{align} & x_{p,d,s} = \begin{cases} 1 & \text{人 } p \text{ が日 } d \text{ のシフト } s \text{ に入る場合} \\ 0 & \text{それ以外} \end{cases} \end{align} $$

制約条件

シフトを作成するための条件から制約条件を定めます。

① 各日の各シフトに必要な人数を満たす制約：

$$ \begin{aligned} \sum_{p \in \text{People}} x_{p,d,s} &= R_s \quad \left( \forall d \in \text{Days}, \forall s \in \text{Shifts} \right) \end{aligned} $$

② 各人の週の勤務回数の最小・最大制約：

\begin{align} \text{Min}_p \leq \sum_{d \in \text{Days}} , \quad \sum_{s \in \text{Shifts}_p} x_{p,d,s} \leq \text{Max}_p \quad \left( \forall p \in \text{People} \right) \end{align}

③ 各人は1日に最大1シフトまで：

$$ \begin{aligned} \sum_{s \in \text{Shifts}_p} x_{p,d,s} \leq 1 \quad \left( \forall p \in \text{People},\ \forall d \in \text{Days} \right) \end{aligned} $$

④ 「遅番→早番」の連続勤務を禁止：

$$ \begin{align} & x_{p,d,\text{late}} + x_{p,d_{\text{next}},\text{early}} \leq 1 \quad \left( \forall p \in \text{People}_{\text{late+early}},\ \forall d \in \text{Days} \right) \\\\ & d_{next} \text{は} d \text{の翌日で、日曜日の翌日は月曜日になります（循環）} \\\\ & \text{People}_{\text{late+early}} \text{：「遅番」と「早番」両方のシフトに入ることができる人たちの集合} \end{align} $$

Pythonで実装する

線形計画問題を実装するために、今回はPythonライブラリのPuLPというソルバーを利用します。
PuLPは厳密解法で解を求めるため、制約条件の変数が多くなりすぎると計算時間が非常に長くなることがありますが、今回はそれほど変数の数が多くないためPuLPで計算をしたいと思います。

ライブラリのインストール

pip install pulp

実際のコード

from pulp import LpProblem, LpVariable, lpSum, LpMinimize, LpBinary, LpStatus, value
import itertools

# 定数の定義
days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
shifts = ['early', 'mid', 'late']
people = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L']

# 各シフトに必要な人数
required = {'early': 2, 'mid': 3, 'late': 2}

# 各人の希望・可能なシフトと勤務日数制約
availability = {
    'A': {'shifts': ['early', 'mid', 'late'], 'min': 5, 'max': 6},
    'B': {'shifts': ['early', 'mid', 'late'], 'min': 4, 'max': 5},
    'C': {'shifts': ['early', 'mid', 'late'], 'min': 3, 'max': 5},
    'D': {'shifts': ['early'], 'min': 3, 'max': 5},
    'E': {'shifts': ['early'], 'min': 2, 'max': 4},
    'F': {'shifts': ['early'], 'min': 2, 'max': 3},
    'G': {'shifts': ['mid', 'late'], 'min': 2, 'max': 3},
    'H': {'shifts': ['mid', 'late'], 'min': 2, 'max': 3},
    'I': {'shifts': ['mid', 'late'], 'min': 4, 'max': 5},
    'J': {'shifts': ['mid', 'late'], 'min': 4, 'max': 5},
    'K': {'shifts': ['late'], 'min': 4, 'max': 5},
    'L': {'shifts': ['late'], 'min': 2, 'max': 6},
}

# 変数の定義
x = LpVariable.dicts("shift", (people, days, shifts), 0, 1, LpBinary)

# 問題の定義（どちらでもよいが、とりあえず最小化問題に設定する）
prob = LpProblem("Shift_Scheduling", LpMinimize)

# 目的関数（ゼロに設定して可行性問題にする）
prob += 0

# 制約条件：各日の各シフトに必要な人数を満たす
for d in days:
    for s in shifts:
        prob += lpSum(x[p][d][s] for p in people if s in availability[p]['shifts']) == required[s]

# 制約条件：各人の週の勤務日数制限
for p in people:
    total_shifts = lpSum(x[p][d][s] for d in days for s in availability[p]['shifts'])
    prob += total_shifts >= availability[p]['min']
    prob += total_shifts <= availability[p]['max']

# 制約条件：各人が1日に1つのシフトしか担当しない
for p in people:
    for d in days:
        prob += lpSum(x[p][d][s] for s in availability[p]['shifts']) <= 1

# 制約条件：「遅番→早番」の禁止
for p in people:
    if 'late' in availability[p]['shifts'] and 'early' in availability[p]['shifts']:
        for i, d1 in enumerate(days):
            d2 = days[(i + 1) % len(days)]  # 次の日（循環）
            prob += x[p][d1]['late'] + x[p][d2]['early'] <= 1

# 解を計算する
prob.solve()

# 結果を整形
schedule = {
    d: {s: [] for s in shifts} for d in days
}

for p in people:
    for d in days:
        for s in shifts:
            # x[p][d][s] に結果が保存されているので、valueで取り出す
            if s in availability[p]['shifts'] and value(x[p][d][s]) == 1:
                schedule[d][s].append(p)

schedule

計算によって得られた解

各人のシフトに入る条件が満たされた状態となっています。

{'Mon': {'early': ['D', 'F'], 'mid': ['A', 'B', 'J'], 'late': ['I', 'K']},
 'Tue': {'early': ['C', 'D'], 'mid': ['A', 'I', 'J'], 'late': ['B', 'L']},
 'Wed': {'early': ['D', 'F'], 'mid': ['A', 'C', 'G'], 'late': ['B', 'K']},
 'Thu': {'early': ['D', 'E'], 'mid': ['A', 'C', 'H'], 'late': ['I', 'K']},
 'Fri': {'early': ['A', 'E'], 'mid': ['B', 'C', 'J'], 'late': ['H', 'I']},
 'Sat': {'early': ['D', 'F'], 'mid': ['A', 'C', 'J'], 'late': ['K', 'L']},
 'Sun': {'early': ['B', 'E'], 'mid': ['G', 'H', 'I'], 'late': ['J', 'K']}}

実装コードの解説

変数の定義

# 変数の定義
x = LpVariable.dicts("shift", (people, days, shifts), 0, 1, LpBinary)

上のコードで使用する変数を定義します。

引数は前から順に、
変数名の prefix（shift）
変数名の postfix のリスト（people, days, shifts）
変数の最小値
変数の最大値
変数の種類（LpBinary：0と1のみ）
となっています。

オブジェクトの定義

# 問題の定義（どちらでもよいが、とりあえず最小化問題に設定する）
prob = LpProblem("Shift_Scheduling", LpMinimize)

上のコードでPuLPのオブジェクトを作成します。

引数は、名前と目的関数の最小化（LpMinimize）または最大化（LpMaximize）を設定します。

目的関数の設定

# 目的関数（ゼロに設定して可行性問題にする）
prob += 0

今回は最小化かつ最大化となるため、目的関数を0に設定します。
ここで注意したいのが、PuLPのオブジェクトに加える ことです。

例えば、x + y を最小化したいときは以下のようにします。

prob = pulp.LpProblem('example', pulp.LpMinimize)
x = pulp.LpVariable('a', 0, 1)
y = pulp.LpVariable('b', 0, 1)
prob += x + y

制約条件の設定

各制約条件を設定しますが、目的関数の設定方法と同様にPuLPのオブジェクトに加えます。

計算の実行

# 解を計算する
prob.solve()

上のコードを実行することで厳密解が得られます。
結果は、設定した変数（x）に保存されており、PuLPライブラリの value メソッドで取り出します。

for p in people:
    for d in days:
        for s in shifts:
            # x[p][d][s] に結果が保存されているので、valueで取り出す
            if s in availability[p]['shifts'] and value(x[p][d][s]) == 1:
                schedule[d][s].append(p)

さいごに

厳密解法では処理時間が長くなることがあるため制約条件を見直すことが必要になるかもしれません。また、近似解法を利用することによって処理時間が短くなることがあります。近似解法については別の記事にできればと思っています。

PuLPで線形計画問題を解くことで、身の回りの問題が解決できそうと感じてもらえたら幸いです。

【参考文献】