Ikegami-2shift-Data1 を正規化して pulp で定式化した話

自分で使うためにやっていたことを，ざっくりとではありますが整理して，[数理最適化Advent Calender 2023](https://qiita.com/advent-calendar/2023/mathematical-optimization)の 6 日目の記事にしました．

# やったこととその背景

最適化のアルゴリズムを開発したりソルバーの動作を確認する際にはベンチマークを使用します．LP であれば [netlib](https://www.netlib.org/lp/data/)，MILP であれば [MIPLIB2017](https://miplib.zib.de/) などが広く知られていますが，これらのベンチマークでは

* 解くことが難しいインスタンスが多く，アルゴリズムの動作確認には適さない
* モデルの意味がわかりづらい

といった問題があり，特にアルゴリズム開発の初期段階では扱いづらいものになっています．

より扱いやすいベンチマークとしては，[池上先生が公開されているナーススケジューリング問題](http://cricket.ikegami-lab.tokyo/~atsuko/DATA/data.html)があり，このモデルは

* 整数線型計画問題(ILP)として定式化可能
* データサイズが小さい
* 簡単すぎない（最適解・実行可能解が自明ではない）
* 難しすぎない[^1]（少なくとも，制約条件を緩和すれば実行可能解を簡単に見つけられる）
* モデルの意味がわかっているので，
    * 自分でモデルを修正して，より簡単な問題・難しい問題を作ることができる
    * アルゴリズムの振る舞いの理由を解釈しやすい
* 実務的なモデルなので解けると嬉しい
* 有名なモデルなので人に説明しやすい

と，アルゴリズムの動作確認や実験に非常に使いやすいものになっています．…が，公開されているデータのフォーマットが ~~今となっては神~~ 昨今のデータサイエンティストのスキルセットとツールでは扱いづらい形式になっているので，[正規化したデータと，そのデータを使ったサンプルモデル](https://github.com/AI448/Normalized_Ikegami-2shift-DATA1)を作りました．

このモデルを LP ファイルや MPS ファイルとして出力する，または（アルゴリズムの開発を python で行うなら）モデルから数式を直接取得するなどして，アルゴリズム開発に使用することができます．

なお，今回作成したデータ・モデルは，作成時のミスや定式化の解釈間違いが発生していないことを十分に確認できておらず，実装していない制約条件もあります．元の論文の結果を再現する目的で使用する場合にはご自身で確認・変更をお願いいたします．

# 作成したデータ

以下のように正規化したテーブルを設計し，各テーブルをシートとして 1 つの Excel ブックにしました．
太字のキーは主キー（テーブルの行を一意に定めるキー）を表します．正規化の程度はいろいろ考えられますが，今回は手修正が容易になることを優先して，一部の小さいデータでは NULL 値を許容してでもテーブル数が少なくなるように設計しています．

## member

メンバー（ナース）に紐づく値を設定

| キー | 意味 |
|:-----------|:------------|
| **member** | メンバー |
| weekends_off_lower | 週末2連休の取得回数下限 |
| weekends_off_upper | 週末2連休の取得回数上限 |

## shift

シフトに紐づく値を設定

| キー | 意味 |
|:-----------|:------------|
| **shift** | シフト種別 |
| consecutive_shifts_lower | shift を連続して割り当てられる日数の下限 |
| consecutive_shifts_upper | shift を連続して割り当てられる日数の上限 |
| days_between_shifts_lower | shift を割り当てられる間隔の下限 |
| days_between_shifts_upper | shift を割り当てられる間隔の上限 |

## ng_pattern

割当不可能な勤務パターンを設定

| キー | 意味 |
|:-----------|:------------|
| **ng_pattern** | パターンを区別するためのID |
| **order** | パターンの何日目の勤務かを表す |
| shift | ng_pattern の order 日目の勤務シフト |

夜勤明け翌日の日勤など，割り当ててはいけないパターンを設定するために使用し，例えば NIGHT（夜勤） -> MORNING(夜勤明けの朝の勤務) -> DAY(日勤)を禁止するのであれば以下のように設定する．

| **ng_pattern** | **order** | shift |
|:-----------|:------------|:----|
| NG1 | 0 | NIGHT |
| NG1 | 1 | MORNING |
| NG1 | 2 | DAY |

## member_shift

メンバーとシフトに紐づく値を設定

| キー | 意味 |
|:-----------|:------------|
| **member** | |
| **shift** | |
| number_of_working_lower | member が shift を行う回数の下限 |
| number_of_working_upper | member が shift を行う回数の上限 |

## group_member
グループとメンバーの対応付を設定

| キー | 意味 |
|:-----------|:------------|
| **member** | |
| **group** | member が属する group |

※ group も主キーである点に注意（=メンバーは複数のグループに属する可能性がある）

## group_date_shift

グループ・日付・シフトに紐づく値を設定

| キー | 意味 |
|:-----------|:------------|
| **group** | |
| **date** | |
| **shift** | |
| number_of_workers_lower | date での， group に属し shift の勤務を行うメンバーの人数の下限 |
| number_of_workers_upper | date での， group に属し shift の勤務を行うメンバーの人数の上限 |

## request_shift

メンバーの希望勤務を設定

| キー | 意味 |
|:-----------|:------------|
| **member** | |
| **date** | |
| shift | member の date での希望シフト |

## refused_shift

メンバー・日付ごとの割当不可シフトを設定

| キー | 意味 |
|:-----------|:------------|
| **member** | |
| **date** | |
| shift | date での member への割当不可シフト |

## past_shift

前月の割当シフト

| キー | 意味 |
|:-----------|:------------|
| **member** | |
| **date** | 前月の日付 |
| shift | member が date に shift の勤務を行ったことを表す |

※今回のサンプルモデルでは使用していません．（...めんどくさくて実装しませんでしたが，実務でのシフトスケジューリングでは，前月末と矛盾しない計画を立てることは重要な要件になります）

# 作成したモデル

定式化の詳細は[こちらのコード](https://github.com/AI448/Normalized_Ikegami-2shift-DATA1/blob/main/sample.py)または挙げている参考文献を見ていただくことにして，ここではデータの読み込みからモデルの構築までの大まかな流れを説明します．

## データの読み込み

```pandas.read_excel``` を使用して，Excel ブックの各シートを ```DataFrame``` として読み込みます．

```python
df_member = pd.read_excel(data_filepath, "member", dtype=str)
```
```python
df_group_member = pd.read_excel(data_filepath, "group_member", dtype=str)
```

## データの構築

```DataFrame``` から最適化に使用するデータを作成します．[^2]

今回のモデルでは，メンバー・グループ・シフトなどの定式化に現れる要素(数式で書いたときに添え字になるもの)を表すための ```dataclass``` を使った型を定義しました．[^3] ただし，シフト種別は既知であるとして，型を ```Enum``` クラスにしています．

```python
@dataclass(frozen=True)
class Group:
    """グループ"""
    id: str
```
```python
@dataclass(frozen=True)
class Member:
    """メンバー"""
    id: str
```
```python
class Shift(Enum):
    """シフト種別"""
    OFF = "OFF"  # 休み
    DAY = "DAY"  # 日勤
    NIGHT = "NIGHT"  # 夜勤
    MORNING = "MORNING"  # 夜勤明けの朝の勤務
    OTHER = "OTHER"  # その他の勤務
```

これらの型を使って ```set```, ```dict``` としてデータを構築します．

```python
    # メンバーの集合を作成
    members = set(map(Member, df_member["member"]))
```
```python
    # グループの集合を作成
    groups = set(map(Group, df_group_member["group"]))
```
```python
    # グループとメンバーの組の集合を作成
    group_members = set(map(
        lambda x: (Group(x[0]), Member(x[1])),
        zip(df_group_member["group"], df_group_member["member"])
    ))
```
```python
    # (グループ, 日付, シフト種別) -> 割当人数下限 の辞書を作成
    number_of_workers_lower = {
        (Group(group), Date(date), Shift(shift)): int(lower)
        for group, date, shift, lower
        in zip(df_group_date_shift["group"], df_group_date_shift["date"], df_group_date_shift["shift"], df_group_date_shift["number_of_workers_lower"])
     }
```

## 定式化

pulp を使って最適化のモデルを作成します．

まずは変数を定義します．シフトスケジューリングでは，ある人に対してある日にシフトを割り当てるか否かを表す 0-1 変数を定義し，割当問題とするのが典型的な定式化方法になります．

```python
    # member に対して date に shift を割り当てるなら 1 そうでなければ 0
    self.x: Dict[Tuple[Member, date, Shift], pulp.LpVariable] = {}
    for member, date, shift in itertools.product(members, dates, Shift):
        self.x[member, date, shift] = pulp.LpVariable(
            name=f"x_{member.id}_{_date_to_id(date)}_{shift.value}",
            cat=pulp.LpBinary
        )
```

変数とさきほど用意したデータとで制約条件を記述します．
ここでは以下の 2 つの制約条件を例にします．
* 各メンバーには同時に 1 つのシフトしか割り当てることができないという自明な制約
* 各日の勤務人数下限制約
    * この制約条件はグループとシフト種別ごとの条件で，「日勤シフトでは合計10人以上」・「夜勤シフトでは合計4人以上」といった全メンバーの合計だけではなく，「夜勤シフトではベテラングループに属する看護師が1人以上」といった条件を設定できるようになっています．

それぞれ以下のように定式化できます．

```python
    # 各 member は各 date においてちょうど 1 つの shift を実行する
    for member, date in itertools.product(members, dates):
        self.problem += (
            pulp.lpSum(self.x[member, date, shift] for shift in Shift) == 1,
            f"c_assignment_{member.id}_{_date_to_id(date)}"
        )
```
```python
    # 勤務人数合計 >= 勤務人数下限    
    for group, date, shift in itertools.product(groups, dates, Shift):
        self.problem += (
            pulp.lpSum(self.x[member, date, shift] for member in group_to_members[group])
            >= number_of_workers_lower[group, date, shift],
            f"c_number_of_workers_lower_{group.id}_{_date_to_id(date)}_{shift.value}"
        )
```

今回の目的は問題を解くことではなく，アルゴリズムの開発に使用するためのデータを作成することなので，モデルをデータとして出力したときにわかりやすいように制約条件にも名前をつけています．pulp を使った通常の開発では，制約条件の名前は必ずしも必要ではありません．

# 得られる LP ファイル

作成したモデルを使って出力した LP ファイルの一部を抜粋しました．

```
c_assignment_10_20230901: x_10_20230901_DAY + x_10_20230901_MORNING
 + x_10_20230901_NIGHT + x_10_20230901_OFF + x_10_20230901_OTHER = 1
```

こちらは上で定義した 1 つ目の制約条件に対応し，「 2023/9/1 に，メンバー 10 に割り当てられるシフトは DAY, MORNING, NIGHT, OFF, OTHER のいずれかである」を表しています．

```
c_number_of_workers_lower_Skilled_20230901_DAY: x_11_20230901_DAY
 + x_12_20230901_DAY + x_1_20230901_DAY + x_20_20230901_DAY
 + x_21_20230901_DAY + x_2_20230901_DAY + x_3_20230901_DAY >= 1
```
また，こちらは 2 つ目の制約条件に対応し，「2023/9/1 に，Skilled グループに属するメンバーで日勤シフトを行う人が 1 人以上」を表しています．

これで，アルゴリズム開発に使えるベンチマーク問題を作成することができました．

# 参考文献

* [データとソース(GitHub)](https://github.com/AI448/Normalized_Ikegami-2shift-DATA1)
* [池上敦子, 丹羽明, 大倉元宏(1996): "我が国におけるナース・スケジューリング問題", オペレー
ションズ・リサーチ, Vol. 41, No. 8, pp. 436-442.](https://orsj.org/wp-content/or-archives50/pdf/bul/Vol.41_08_436.pdf)
* http://cricket.ikegami-lab.tokyo/~atsuko/DATA/data.html


[^1]: 2023年現在の計算機と最適化ソルバーであれば解くことができますが，発表当時は難しかったものと思われます．
[^2]: ```DataFrame``` を使って最適化モデルを記述することも可能ですが，複雑になりがちなので，```dict```, ```set``` などの python の組み込み型を使う方をおすすめします．
[^3]: 型アノテーションはシステム開発では有用ですが，今回は開発のためのベンチマーク問題の作成が目的なので，凝りすぎているかもしれません．