遺伝的アルゴリズムでナーススケジューリング問題（シフト最適化）を解く

この記事は
===

* ナーススケジューリング問題という最適化の問題を遺伝的アルゴリズムで解いてみたらまあまあの精度が出たので記録です
* pythonのdeapというライブラリを使っています

前提
===

ナーススケジューリング問題
---

* [ナーススケジューリング問題](http://www.msi.co.jp/nuopt/glossary/term_9ad8057bc9aa1325cae9843e345fd8c9034ef119.html)というものがあります

> 病院等の医療施設に勤める看護師の勤務スケジュールを決定する問題のことであり，シフトスケジューリング問題の代表例である．日勤・夕勤・夜勤等の複雑なシフト勤務や多岐に渡る制約の考慮のため，実際にスケジュールを求めるのは人手では手間のかかる困難な作業であり・・・

* 要するにシフト勤務のスケジュールを自動的に最適に組むというものです
* 制約が色々あって、必要人数を満たすような基本的なものから、公平性、必須な資格/役割、この２人は（仲が悪いから？）一緒に入れないなど、いくらでも複雑になりうる問題て感じで、完璧な解答を作るのは困難なので、近似解を求めるような分野みたいです
*

遺伝的アルゴリズム
---

* 上記問題を解くやり方はいろいろあるみたいなのですが、[遺伝的アルゴリズム](https://ja.wikipedia.org/wiki/%E9%81%BA%E4%BC%9D%E7%9A%84%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0)を使う方法が割とメジャーらしいです
* 遺伝的アルゴリズムは大体下記のような手法で最適化を行うものです
  * データ（遺伝子らしく通常は1次元配列orタプル）を内包する「個体」をたくさんつくる
      * 初めに作る個体のデータはランダム
  * 2つの個体を組み合わせる（交叉）、ランダムに変異する（突然変異）、などで個体に変化を与えながら、世代を進行させていく
  * 予め定義されている制約を個体のデータがどれだけ満たしているかどうかで個体のスコアをつける
      * 優秀なスコアの個体ほど世代進行しても生き残るように調整する
  * 一定の世代を経過させ、全体で最もスコアの高かった個体を近似解として得る
* 最適化の中でも離散的な解を求める、組み合わせ最適化全般に使える手法です
* 詳細は[こちら](http://algoful.com/Archive/Algorithm/ga)だとかがわかりやすいです
  * 特に1が多いほどハイスコアという制約の「OneMax問題」というのが一番シンプルで分かりやすいです

で、どう適用するのか
===

* 上記の遺伝的アルゴリズムをどう適用するとナーススケジューリング問題が解けるのか？ですが、個体に持たせるデータをシフトのコマにするイメージです。
* 例）たとえば3人のシフトを下記のように組んでみたとします（1がアサイン対象のコマとする）

曜日|月|月|月|火|火|火|水|水|水|木|木|木|金|金|金|土|土|土|日|日|日
:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:
時間帯|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜
従業員0|0|0|0|1|0|0|1|0|0|0|0|0|1|0|0|0|0|0|1|0|0
従業員1|1|1|0|0|0|0|0|1|1|0|0|0|0|0|1|0|0|0|0|0|0
従業員2|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1|1|0|1|0

* ここで従業員0から従業員2までのゼロ・イチを1つの配列につなげていったものを個体のデータとします
* つまり上記の場合、下記の1次元配列データを持たせるわけです

```
[0,0,0,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,
 1,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,0,
 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,0]
```

* この個体が各種制約をどれだけ満たしているかを計算するスコアリングのロジックを適切に組み、遺伝的アルゴリズムのルールに則って世代を進めていきます
* そうすると、より優れた個体が残るというわけです（=より最適な組み合わせが得られる）

Python/deapで実装
===

具体的に実装を進めてみます。

deapの利用
---

* 実装は機械学習関連に強いpythonを利用しました
* [deap](https://github.com/DEAP/deap)という遺伝的アルゴリズムのライブラリが充実しています
* 普通にpipでインストールします（詳細は割愛）

```
pip install deap
```

* 使い方は[こちら](http://testpy.hatenablog.com/entry/2017/01/09/213846)のOneMax問題の例がわかりやすいです

サンプルの仕様
===

今回のサンプルの仕様は下記とします。

従業員とシフト希望
---

* 従業員は10人とし、下記のシフト希望を出していることとします

曜日|月|月|月|火|火|火|水|水|水|木|木|木|金|金|金|土|土|土|日|日|日
:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:
時間帯|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜
必要人数|2|3|3|2|3|3|2|3|3|1|2|2|2|3|3|2|4|4|2|4|4
従業員0|○| | |○| | |○| | |○| | |○| | |○| | |○| | 
従業員1|○|○|○| | | |○|○|○| | | |○|○|○| | | | | | 
従業員2| | | | | | | | | | | | | | | |○|○|○|○|○|○
従業員3|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○
従業員4| | |○| | |○| | |○| | |○| | |○| | |○| | |○
従業員5|○|○|○|○|○|○|○|○|○|○|○|○|○|○|○| | | | | | 
従業員6| | | | | | | | | | | | |○|○|○|○|○|○|○|○|○
従業員7| |○| | |○| | |○| | |○| | |○| | |○| | |○| 
従業員8| | |○| | |○| | |○| | |○| | |○| | |○| | |○
従業員9| | | | | | | | | |○|○|○|○|○|○|○|○|○|○|○|○

* 必要人数というのは、このコマに必要な従業員の人数です

制約
---

下記の制約を満たすほどにスコアが上がるようにします。括弧内の数字は重みづけで、制約の優先度を表します。

* 必要人数とアサイン従業員人数が一致していること(10)
* 従業員が希望していないコマへのアサインはしないこと(100)
* 従業員各々が希望したコマ数の最低1/2はアサインすること(1)
  * せっかく希望してくれたのでできるだけ入れてあげたいという意図
* 管理者は必ず1人配置すること(100)
  * 管理者は従業員3, 5, 9とします
* 同一従業員は1日2コマまでのアサインとすること(10)


コード概要
===

* プログラム全体は[こちら](https://github.com/shouta-dev/nurse-scheduling-ga/blob/master/nurse_scheduling_by_ga.py)
* かなり[こちらのOneMax問題](http://tech.mof-mof.co.jp/blog/ga-one-max-problem.html)のコードを参考（≒コピペ）させていただきました

ポイント
---

* 上記のOneMax問題からナーススケジューリング問題向けに改造したポイントだけ説明しておきます

### 従業員の定義

* まずオブジェクト指向っぽく従業員クラスを定義してます

```python
class Employee(object):
  def __init__(self, no, name, age, manager, wills):
    self.no = no
    self.name = name
    self.age = age
    self.manager = manager

    # willは曜日_時間帯。1は朝、2は昼、3は夜。
    # 例）mon_1は月曜日の朝
    self.wills = wills
...
```

* インスタンス作ってるところは下記

```python
# 朝だけ
e0 = Employee(0, "山田", 40, False, ['mon_1', 'tue_1', 'wed_1', 'thu_1', 'fri_1', 'sat_1', 'sun_1'])

# 月・水・金
e1 = Employee(1, "鈴木", 21, False, ['mon_1', 'mon_2', 'mon_3', 'wed_1', 'wed_2', 'wed_3','fri_1', 'fri_2', 'fri_3'])

# 週末だけ
e2 = Employee(2, "佐藤", 18, False, ['sat_1', 'sat_2', 'sat_3', 'sun_1', 'sun_2', 'sun_3'])

# どこでもOK
e3 = Employee(3, "田中", 35, True, ['mon_1', 'mon_2', 'mon_3',
                                     'tue_1', 'tue_2', 'tue_3',
                                     'wed_1', 'wed_2', 'wed_3',
                                     'thu_1', 'thu_2', 'thu_3',
                                     'fri_1', 'fri_2', 'fri_3',
                                     'sat_1', 'sat_2', 'sat_3',
                                     'sun_1', 'sun_2', 'sun_3'
                                    ])

...

employees = [e0, e1, e2, e3, e4, e5, e6, e7, e8, e9]
```

### シフト（コマ）の定義

* シフトというかコマの集合のクラスを定義しています
* このクラスの「list」というタプルこそが遺伝的アルゴリズムの「個体」にあたります

```python
# シフトを表すクラス
# 内部的には 3(朝昼晩) * 7日 * 10人 = 210次元のタプルで構成される
class Shift(object):
  # コマの定義
  SHIFT_BOXES = [
    'mon_1', 'mon_2', 'mon_3',
    'tue_1', 'tue_2', 'tue_3',
    'wed_1', 'wed_2', 'wed_3',
    'thu_1', 'thu_2', 'thu_3',
    'fri_1', 'fri_2', 'fri_3',
    'sat_1', 'sat_2', 'sat_3',
    'sun_1', 'sun_2', 'sun_3']

  # 各コマの想定人数
  NEED_PEOPLE = [
    2,3,3,
    2,3,3,
    2,3,3,
    1,2,2,
    2,3,3,
    2,4,4,
    2,4,4]

  def __init__(self, list):
    if list == None:
      self.make_sample()
    else:
      self.list = list
    self.employees = []

  # ランダムなデータを生成
  def make_sample(self):
    sample_list = []
    for num in range(210):
      sample_list.append(random.randint(0, 1))
    self.list = tuple(sample_list)
...
```

* 上記の通り、最初はランダムなデータを持つ個体を作っています
* Shiftクラスはもっと長いのですが、後は延々と条件にあったコマとかユーザを取得するデータマニピュレーション的なメソッドが続いているだけです

### スコアリングと重みづけ

* deapの機能を使って、個体のスコアリングを重みづけ付きで実施しています

```python
creator.create("FitnessPeopleCount", base.Fitness, weights=(-10.0, -100.0, -1.0, -100.0, -10.0))
creator.create("Individual", list, fitness=creator.FitnessPeopleCount)

...

def evalShift(individual):
  s = Shift(individual)
  s.employees = employees

  # 想定人数とアサイン人数の差
  people_count_sub_sum = sum(s.abs_people_between_need_and_actual()) / 210.0
  # 応募していない時間帯へのアサイン数
  not_applicated_count = s.not_applicated_assign() / 210.0
  # アサイン数が応募数の半分以下の従業員数
  few_work_user = len(s.few_work_user()) / 10.0
  # 管理者が１人もいないコマ数
  no_manager_box = len(s.no_manager_box()) / 21.0
  # 朝・昼・夜の全部にアサインされている
  three_box_per_day = len(s.three_box_per_day()) / 70.0
  return (not_applicated_count, people_count_sub_sum, few_work_user, no_manager_box, three_box_per_day)

toolbox.register("evaluate", evalShift)
```

* ```return (not_applicated_count, people_count_sub_sum, few_work_user, no_manager_box, three_box_per_day)```に対して```weights=(-10.0, -100.0, -1.0, -100.0, -10.0)```という重みをつけています（順序で指定する）
* ちなみにスコアは小さいほどよいという仕様のようなので、重みはマイナス値にする必要があります

### 進化の回数等を指定

* 下記の通り、300個体で約500世代まで計算させてみます

```python
    pop = toolbox.population(n=300)
    CXPB, MUTPB, NGEN = 0.6, 0.5, 500 # 交差確率、突然変異確率、進化計算のループ回数
```

実行
===

* 実行すると下記のように学習の過程が出力されます

```
$ python nurse_scheduling_by_ga.py 
進化開始
  300 の個体を評価
-- 0 世代 --
  245 の個体を評価
* パラメータ1
  Min 0.242857142857
  Max 0.242857142857
  Avg 0.242857142857
  Std 2.2660056242e-08
* パラメータ2
  Min 0.247619047619
  Max 0.247619047619
  Avg 0.247619047619
  Std 9.49766396283e-09
* パラメータ3
  Min 0.0
  Max 0.0
  Avg 0.0
  Std 0.0
* パラメータ4
  Min 0.142857142857
  Max 0.142857142857
  Avg 0.142857142857
  Std 1.66600046866e-08
* パラメータ5
  Min 0.114285714286
  Max 0.114285714286
  Avg 0.114285714286
  Std 1.19267483008e-08
-- 1 世代 --
  235 の個体を評価
* パラメータ1
  Min 0.238095238095
  Max 0.238095238095
  Avg 0.238095238095
  Std 2.45699769971e-08
* パラメータ2
  Min 0.228571428571
  Max 0.228571428571
  Avg 0.228571428571
  Std 2.38534966016e-08
* パラメータ3
  Min 0.1
  Max 0.1
  Avg 0.1
  Std 1.31048702444e-08
* パラメータ4
  Min 0.047619047619
  Max 0.047619047619
  Avg 0.047619047619
  Std 4.46646616171e-09
* パラメータ5
  Min 0.1
  Max 0.1
  Avg 0.1
  Std 1.31048702444e-08

(snip)

-- 499 世代 --
  219 の個体を評価
* パラメータ1
  Min 0.0333333333333
  Max 0.0333333333333
  Avg 0.0333333333333
  Std 2.98168707519e-09
* パラメータ2
  Min 0.0714285714286
  Max 0.0714285714286
  Avg 0.0714285714286
  Std 8.33000234328e-09
* パラメータ3
  Min 0.1
  Max 0.1
  Avg 0.1
  Std 1.31048702444e-08
* パラメータ4
  Min 0.0952380952381
  Max 0.0952380952381
  Avg 0.0952380952381
  Std 8.93293232343e-09
* パラメータ5
  Min 0.0428571428571
  Max 0.0428571428571
  Avg 0.0428571428571
  Std 4.56253018749e-09
-- 進化終了 --
最も優れていた個体: [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1], (0.0, 0.0, 0.2, 0.047619047619047616, 0.1)
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0
0,1,0,0,0,0,1,1,0,0,0,0,1,1,1,0,0,0,0,0,0
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1
1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,1,1,0,0,1
0,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0
1,1,1,1,1,1,1,1,1,0,0,1,0,1,0,0,0,0,0,0,0
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,1,1
0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,1,0,0,1,0
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,1,0,1,1
0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	1	0	0
0	1	0	0	0	0	1	1	0	0	0	0	1	1	1	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1
1	0	1	0	1	1	0	0	1	0	1	0	1	0	1	0	1	1
0	0	0	0	0	1	0	0	1	0	0	1	0	0	1	0	0	0
1	1	1	1	1	1	1	1	1	0	0	1	0	1	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	1	1
0	1	0	0	1	0	0	1	0	0	1	0	0	0	0	0	1	0
0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	1	1	1
```

* 結果が配列形式、CSV形式、TSV形式で標準出力される手抜き仕様です
* 適当に加工してExcelに張りつけて利用します
* ちなみに実行時間は2分弱位でした

結果
===

* で、なんどか実施してみて、結果は下記となりました
* （上記の実施時の結果とはデータが別物ですが、気にしないでください）

**曜日**|**月**|**月**|**月**|**火**|**火**|**火**|**水**|**水**|**水**|**木**|**木**|**木**|**金**|**金**|**金**|**土**|**土**|**土**|**日**|**日**|**日**|**アサイン数**|**希望数**|**アサイン率**
:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:|:-----:
時間帯|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜|朝|昼|夜| | | 
必要人数|2|3|3|2|3|3|2|3|3|1|2|2|2|3|3|2|4|4|2|4|4| | | 
承認人数|2|3|3|2|3|3|2|3|3|1|2|2|2|3|3|2|4|4|2|4|4| | | 
差分|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0| | | 
 | | | | | | | | | | | | | | | | | | | | | | | | 
従業員0|1|0|0|1|0|0|1|0|0|0|0|0|0|0|0|1|0|0|1|0|0|5|7|71.4%
従業員1|0|1|1|0|0|0|0|1|1|0|0|0|0|1|1|0|0|0|0|0|0|6|9|66.7%
従業員2|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1|1|0|1|1|4|6|66.7%
[管]従業員3|0|1|0|0|1|0|0|1|0|1|1|0|1|0|0|0|0|0|0|1|0|7|21|**33.3%**
従業員4|0|0|1|0|0|1|0|0|1|0|0|0|0|0|0|0|0|1|0|0|1|5|7|71.4%
[管]従業員5|1|0|**1**|**1**|**1**|1|1|0|0|0|0|1|0|0|1|0|0|0|0|0|0|8|15|53.3%
従業員6|0|0|0|0|0|0|0|0|0|0|0|0|1|1|0|1|1|1|0|0|0|5|9|55.6%
従業員7|0|1|0|0|1|0|0|1|0|0|1|0|0|0|0|0|1|0|0|1|0|6|7|85.7%
従業員8|0|0|0|0|0|1|0|0|1|0|0|0|0|0|1|0|0|0|0|0|1|4|7|57.1%
[管]従業員9|0|0|0|0|0|0|0|0|0|0|0|1|0|1|0|0|1|1|**1**|**1**|**1**|7|12|58.3%
 | | | | | | | | | | | | | | | | | | | | | | | | 
管理者数|1|1|1|1|2|1|1|1|**0**|1|1|2|1|1|1|**0**|1|1|1|2|1| | | 

* 太字が条件をみたせなかったところです
* かなり複雑な条件を指定した割りには大体うまくいってます

終わりに
===

* ここまでできてれば、人手で微調整すればすぐにシフト作れそうです
* ただ、2分弱は時間かかりすぎだなーと思いました。
  * [scoop](https://github.com/soravux/scoop)で分散化したりしていくつか高速化を試みたのですが結局CPUネックになってしまってあまり早くなりませんでした
  * DeepLearningみたいにGPU使ったら早くなるかもしれません
      * [論文レベル](http://mikilab.doshisha.ac.jp/dia/monthly/monthly2011/mlm126/ykurano/ykurano.pdf)では出てきましたが、実装は見つかりませんでした。。cudaとか使って自作するのはきつい。。
* ちなみに、私ごとですが姉が看護師で、シフト作るのに一週間くらいスケジュールとにらめっこしているらしいので、本気で制約条件作り込めば結構役に立つんじゃないのかなーと思いました