航空機における搭乗戦略の最適化 - OR機関誌10月号より

# これなに
[OR学会](http://www.orsj.or.jp/)(問題解決学である[オペレーションズ・リサーチ](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%83%9A%E3%83%AC%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3%E3%82%BA%E3%83%BB%E3%83%AA%E3%82%B5%E3%83%BC%E3%83%81)の研究者の集まり)の機関誌10月号は、[「**学生たちのOR**」特集](http://www.orsj.or.jp/e-library/elcorsj.html#6110)となっており、大学生が取組んだ様々な卒業論文や修士論文などの要旨が30編、紹介されています。

この中から、適当に最適化の問題をピックアップしてPythonで解いてみたいと思います。
準備として、pandas, pulp, ortoolpy が必要です。環境構築は、[組合せ最適化を使おう](http://qiita.com/Tsutomu-KKE@github/items/bfbf4c185ed7004b5721#%E3%82%BD%E3%83%95%E3%83%88%E3%81%AE%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB)を参考にしてください。

# 航空機における搭乗戦略の最適化
論文「[航空機における搭乗戦略の最適化](http://www.orsj.or.jp/archive2/or61-10/or61_10_686.pdf)」の問題を使わせてもらいましょう。

> 乗客6人を3グループに分ける。グループ1から順番に搭乗する。同一グループ内は、ランダムな順番とする。総搭乗時間を最小化したい。

# 考え方

座席$i$の乗客が先に搭乗し、座席$j$の乗客が後に乗った場合の混雑度$a_{ij}$を要素とする行列$A$が与えられるものとします。ここでは、単純に、この混雑度の総和を最小化(=混雑しない組合せの混雑度の最大化)をしましょう。また、同一グループでは、混雑度を半分とします。

混雑しない組合せの混雑度 = (j,iの順番に搭乗する場合の$a_{ij}$の和) + (i,jが同一グループの場合の$a_{ij}$の和)/2

# Pythonで解く
まず、ランダムなデータを作成します。

```py3:python
import numpy as np, pandas as pd
from pulp import *
from ortoolpy import addvar, addvars

m, n = 3, 6 # グループ数と乗客数
rm, rn = range(m), range(n)
np.random.seed(1)
A = np.random.rand(n, n).round(3)
A[np.diag_indices(n)] = 0
A
>>>
array([[ 0.   ,  0.72 ,  0.   ,  0.302,  0.147,  0.092],
       [ 0.186,  0.   ,  0.397,  0.539,  0.419,  0.685],
       [ 0.204,  0.878,  0.   ,  0.67 ,  0.417,  0.559],
       [ 0.14 ,  0.198,  0.801,  0.   ,  0.313,  0.692],
       [ 0.876,  0.895,  0.085,  0.039,  0.   ,  0.878],
       [ 0.098,  0.421,  0.958,  0.533,  0.692,  0.   ]])
```

座席(Pos)の乗客がグループ(Group)かどうかを管理する変数表を作成します。

```py3:python
tg = pd.DataFrame(((i, j+1) for i in rn for j in rm), columns=['Pos', 'Group'])
tg['Var'] = addvars(len(tg), cat=LpBinary)
tg[:3]
```
<table>
  <thead>
      <th></th>
      <th>Pos</th>
      <th>Group</th>
      <th>Var</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>1</td>
      <td>v1</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>2</td>
      <td>v2</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0</td>
      <td>3</td>
      <td>v3</td>
    </tr>
  </tbody>
</table>

混雑度Aがかかるかどうかを管理する変数表を作成します。
VarNは、座席Firstに後に搭乗し、座席Secondに先に乗った場合を、
VarHは、座席Firstと座席Secondが同一グループの場合を表すものとします。

```py3:python
tp = pd.DataFrame(((i, j, A[i,j]) for i in rn for j in rn if A[i,j]),
    columns=['First', 'Second', 'A'])
tp['VarN'] = addvars(len(tp), cat=LpBinary)
tp['VarH'] = addvars(len(tp), cat=LpBinary)
tp[:3]
```

<table>
  <thead>
    <tr>
      <th></th>
      <th>First</th>
      <th>Second</th>
      <th>A</th>
      <th>VarN</th>
      <th>VarH</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>1</td>
      <td>0.720</td>
      <td>v19</td>
      <td>v48</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>3</td>
      <td>0.302</td>
      <td>v20</td>
      <td>v49</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0</td>
      <td>4</td>
      <td>0.147</td>
      <td>v21</td>
      <td>v50</td>
    </tr>
  </tbody>
</table>

定式化して解きます。制約式の説明は、面倒なので省略。

```py3:python
m = LpProblem(sense=LpMaximize) # 数理モデル
m += lpDot(tp.A, tp.VarN) + lpDot(tp.A, tp.VarH) / 2 # 目的関数
for i in rn:
    m += lpSum(tg[tg.Pos == i].Var) == 1 # 必ず何れかのグループに所属する
for _, r in tp.iterrows():
    tf = tg[tg.Pos == r.First]
    ts = tg[tg.Pos == r.Second]
    m += (lpDot(tf.Group, tf.Var) - lpDot(ts.Group, ts.Var) - 1)/n + 1 >= r.VarN
    m += (lpDot(tf.Group, tf.Var) - lpDot(ts.Group, ts.Var))/(n-1) + 1 >= r.VarH
    m += (lpDot(ts.Group, ts.Var) - lpDot(tf.Group, tf.Var))/(n-1) + 1 >= r.VarH
m.solve() # ソルバー(CBC)の実行
tg['Val'] = tg.Var.apply(value) # 結果
tg[tg.Val > 0] # 解の表示
```

<table>
  <thead>
    <tr>
      <th></th>
      <th>Pos</th>
      <th>Group</th>
      <th>Var</th>
      <th>Val</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>1</th>
      <td>0</td>
      <td>2</td>
      <td>v2</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>1</td>
      <td>2</td>
      <td>v5</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>8</th>
      <td>2</td>
      <td>3</td>
      <td>v9</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>9</th>
      <td>3</td>
      <td>1</td>
      <td>v10</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>14</th>
      <td>4</td>
      <td>3</td>
      <td>v15</td>
      <td>1.0</td>
    </tr>
    <tr>
      <th>15</th>
      <td>5</td>
      <td>1</td>
      <td>v16</td>
      <td>1.0</td>
    </tr>
  </tbody>
</table>

座席(Pos)ごとのグループ(Group)が求まりました。
このアプローチは、離散変数が多いので、乗客数が増えると、計算時間が爆発します。その場合は、元の論文のように局所探索法などの近似解法が有効でしょう。

# ORセミナーの紹介
11/12(土)に開かれる[第3回ORセミナー](http://www.orsj.or.jp/activity/seminar.html#semi3)(Python言語によるビジネスアナリティクス)では、最適化・統計分析・機械学習などのオペレーションズ・リサーチの分野で必要なツールの紹介があります。
このセミナーの参加者特典として、2016年度と2017年度の年会費＋入会費が免除で、学会正会員になれます。正会員になると、上記、機関誌(年12冊)や論文誌を受取れたり、シンポジウムの無料参加などの特典があります。

以上