Python mipライブラリを使った数理最適化 応用編

## はじめに

　前回の記事[Python mipライブラリを使った数理最適化](https://qiita.com/miso_taku/items/be16af812f3d81486bd4)では数理最適化の基礎を紹介しました。本記事ではその応用例を取り上げます。


## 原材料輸送費の最適化

### 問題背景
　プリン専門店チェーンは、効率的な物流を目指し、砂糖・卵・牛乳といった原材料の配送コストを最小限に抑えたいと考えています。それぞれの専門店は異なる需要量を抱え、原材料の供給元である工場（養鶏場）も生産量に制限があります。また、配送にはコストがかかり、各工場から各店舗への配送コストも異なります。このような状況で、チェーン全体としての物流コストを最小化する方法を見つけ出します。

### 問題の概要
チェーン店舗とその需要量
各店舗のプリン需要量は以下の通りです：

| 店舗 | 需要量 | 
| - | - |
| 店A | 500 |
| 店B | 400 |
| 店C | 300 |

原材料の工場と生産量
プリンを作るために必要な原材料は「卵」「牛乳」「砂糖」の3種類です。各工場の生産量は次の通りです：

- 養鶏場（卵）

| 工場 | 生産量 | 
| - | - |
| 養鶏場a | 700 | 
| 養鶏場b | 500 |


- 牛乳工場

| 工場 | 生産量 | 
| - | - |
| 牛乳工場c | 900 | 
| 牛乳工場d | 300 | 


- 砂糖工場

| 工場 | 生産量 | 
| - | - |
| 砂糖工場e | 800 | 
| 砂糖工場f | 400 |

必要な原材料
プリンを1作るには以下の原材料が卵、牛乳、砂糖が1ずつ必要になります: 

配送コスト
各工場から各店舗への配送コスト（1あたり）は次の表のように定義されています：

| 工場 | 店A | 店B | 店C |
| - | - | - | - |
| 養鶏場a | 2 | 2 | 8 |
| 養鶏場b | 4 | 4 | 4 | 
| 牛乳工場c | 5 | 5 | 2 |
| 牛乳工場d | 5 | 4 | 3 |
| 砂糖工場e	| 3	| 4	| 5 |
| 砂糖工場f	| 6	| 3	| 3 | 

### 目的

　各店舗の需要を満たすために必要な原材料を工場から配送する際、総配送コストを最小化することを目指します。

### コード

　今回もgoogle colabで実施するので、セルで下記の内容を実施して下さい。

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

```python
!pip install mip
```

#### ライブラリのインポート

```python
from mip import Model, xsum, minimize, CONTINUOUS
```


#### 各種変数の定義

```python
# 店舗のリストと需要量
demand = {"A": 500, "B": 400, "C": 300}

# 各工場の生産量
supply = {
    "egg": {"a": 700, "b": 500},
    "milk": {"c": 900, "d": 300},
    "sugar": {"e": 800, "f": 400}
}

# 配送コスト
cost = {
    'egg': {
        'a': {'A': 2, 'B': 2, 'C': 8}, 
        'b': {'A': 4, 'B': 4, 'C': 4}
        },
    'milk': {
        'c': {'A': 5, 'B': 5, 'C': 2}, 
        'd': {'A': 5, 'B': 4, 'C': 3}
        },
    'sugar': {
        'e': {'A': 3, 'B': 4, 'C': 5}, 
        'f': {'A': 6, 'B': 3, 'C': 3}
        }
}

```

#### モデルの作成と変数の設定

```python
# 最適化モデルの作成
model = Model()

# 変数の作成
# x:各店舗への配送料を表す変数
x = {}
for material in supply:
    x[material] = {}
    for factory in supply[material]:
        x[material][factory] = {}
        for store in demand:
            x[material][factory, store] = model.add_var(var_type=CONTINUOUS, name=f"{material}_{factory}_{store}", lb=0)
```

#### 目的関数の設定

```python
model.objective = minimize(
    xsum(
        cost[material][factory][store] * x[material][factory, store]
        for material in supply
        for factory in supply[material]
        for store in demand
    )
)
```

#### 制約条件の追加

```python
# 店舗ごとに需要量を満たす
for store in demand:
    for material in supply:
        model += xsum(x[material][factory, store] for factory in supply[material]) >= demand[store]

# 工場の生産量の制約
for material in supply:
    for factory in supply[material]:
        model += xsum(x[material][factory, store] for store in demand) <= supply[material][factory]
```

#### 最適化の実行と結果表示

```python
# 最適化の実行
model.optimize()

# 結果の表示
if model.num_solutions:
    for material in supply:  # materialでループ
        for factory in supply[material]:  # factoryでループ
            for store in demand:  # storeでループ
                quantity = x[material][factory, store].x  # タプルキーを使用して変数にアクセス
                if quantity > 0:
                    if material == 'egg':
                        # 養鶏場からの配送を表示
                        print(f"養鶏場 {factory} から店舗 {store}: {quantity}")
                    elif material == 'milk':
                        # 牛乳工場からの配送を表示
                        print(f"牛乳工場 {factory} から店舗 {store}: {quantity}")
                    elif material == 'sugar':
                        # 砂糖工場からの配送を表示
                        print(f"砂糖工場 {factory} から店舗 {store}: {quantity}")
# トータルの配送コストを表示
print(f"総配送コスト: {model.objective_value}")

# ※コードのネストが深くて、あんまり良くないコードやな…
```


出力

```出力
養鶏場 a から店舗 A: 300.0
養鶏場 a から店舗 B: 400.0
養鶏場 b から店舗 A: 200.0
養鶏場 b から店舗 C: 300.0
牛乳工場 c から店舗 A: 500.0
牛乳工場 c から店舗 B: 100.0
牛乳工場 c から店舗 C: 300.0
牛乳工場 d から店舗 B: 300.0
砂糖工場 e から店舗 A: 500.0
砂糖工場 e から店舗 B: 300.0
砂糖工場 f から店舗 B: 100.0
砂糖工場 f から店舗 C: 300.0
総配送コスト: 12100.0
```

## スタッフのシフトスケジュールの最適化

### 問題設定
　シフト作成は、スタッフの希望を最大限考慮しつつ、職場の運営を効率化するための重要なタスクです。しかし、希望を満たしながら、全体のバランスを取るのは容易ではありません。そこで、Pythonのmipライブラリを用いて、希望シフトの成立数を最大化するシフト作成問題を解く方法を解説します。

### 目標
　今回の目標は、プリン専門店に務めるスタッフが希望するシフトをできるだけ多く成立させつつ、以下の業務要件を満たすシフト表を作成することです。

### スタッフの人数
スタッフは全員で5人（伊藤, 黒澤, 成瀬, 響野, 久遠）です。

### シフトの種類と曜日
シフトには以下の2種類があります：
- 朝番
- 夜番

対象期間は1週間（月曜日～日曜日）とします。

### 制約条件
シフトを割り当てる際には、以下の制約を満たす必要があります：

- 1日1シフト制約
    - 各スタッフは、1日のうち1つのシフト（朝番または夜番）のみ担当します
- 休日日数制約
    - 各スタッフは、1週間のうち2～3日は休みを取らなければなりません
- 必要人数制約
    - 各曜日のシフトには以下の人数が必要です：
        - 朝番：1人
        - 夜番：2人


### コード　
　今回もgoogle colabで実施するので、セルで下記の内容を実施して下さい。

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

```python
!pip install mip
```

#### ライブラリのインポート

```python
from mip import Model, maximize, xsum
```

#### 各種変数の定義

```python
staffs = ["伊藤", "黒澤", "成瀬", "響野", "久遠"]
days = ["月", "火", "水", "木", "金", "土", "日"]
shifts = ["朝", "夜", "休"]
n0 = len(staffs)
n1 = len(days)
n2 = len(shifts)
```

#### スタッフの希望休日設定

```python
days_off_desired_staff = [
    ["土", "日"],      # 伊藤の希望休
    ["月", "水", "金"], # 黒澤の希望休
    ["土", "日"],      # 成瀬の希望休
    ["月", "火"],      # 響野の希望休
    ["水", "木", "金"] # 久遠の希望休
]

day_to_index = {day: i for i, day in enumerate(days)}
wish_days = [[day_to_index[day] for day in staff_wish] for staff_wish in days_off_desired_staff]
```

#### モデルの作成と変数の定義

```python
model = Model()
x = model.add_var_tensor((n0, n1, n2), "x", var_type="B")
```

#### 目的関数の設定

```python
x_wish_days = [] # 休みの希望の変数のリスト
for x_i, wish_day in zip(x, wish_days):
    for day in wish_day:
        # 該当スタッフの曜日ごとの休み変数を追加
        x_wish_days.append(x_i[day, 2])

model.objective = maximize(xsum(x_wish_days))
```

#### 制約条件の定義

```python
# スタッフごと
for i0 in range(n0):
    # 日番号ごと
    for i1 in range(n1):
        # 1つのシフトを割り当て
        model += xsum(x[i0][i1]) == 1
    model += xsum(x[i0, :, 2]) >= 2 # 休み2日以上
    model += xsum(x[i0, :, 2]) <= 3 # 休み3日以下

# 曜日ごと
for i1 in range(n1):
    model += xsum(x[:, i1, 0]) == 1 # 朝1人
    model += xsum(x[:, i1, 1]) == 2 # 夜2人
```

#### 最適化の実行

```python
model.verbose = 0
model.optimize()
```

#### 結果の表示

```python
# 結果表示関数
def result_display(result):
    print("     月火水木金土日")
    for i0 in range(n0):
        print(f"{staffs[i0]}:", end="")
        for i1 in range(n1):
            print(shifts[result[i0, i1]], end="")
        print()

# 結果表示
if model.status.value == 0:
    val = x.astype(float, subok=False).round().astype(int)
    result = ([0, 1, 2] * val).sum(axis=2)
    result_display(result)
else:
    print("最適化失敗")
```

結果出力
```
     月火水木金土日
伊藤:朝休朝夜朝休休
黒澤:休夜休夜休朝朝
成瀬:夜朝夜休夜休休
響野:休休夜朝夜夜夜
久遠:夜夜休休休夜夜
```

## まとめ

　この記事では、Pythonのmipライブラリを用いて、輸送費の最適化問題とスタッフの希望を考慮したシフト作成問題を解決しました。このように数理最適化を使えばいろいろな問題に対して最適な解を見つけることができます。

 ## おまけ

 mipライブラリのよく使う関数や処理など一覧

| 関数や処理 | 説明 |
| --- | --- |
| model = Model() | モデルの作成 |
| model.add_var(変数) | 変数を追加 |
| model.add_var_tensor((変数))、変数 | 多次元配列の形で変数を作成します |
| model.objective　= (成立させる計算式)| 目的関数の設定または取得 |
| model += xsum(変数ベクトル) | 制約条件の追加 |
| model.optimize() | 最適化を実行 |
| model.status | 最適化ステータス(実行結果)を取得 |
| model.objective_value | 最適化後の目的関数の値を取得 |