これなに
OR学会(問題解決学であるオペレーションズ・リサーチの研究者の集まり)の機関誌10月号は、「学生たちのOR」特集となっており、大学生が取組んだ様々な卒業論文や修士論文などの要旨が30編、紹介されています。
この中から、適当に最適化の問題をピックアップしてPythonで解いてみたいと思います。
準備として、pandas, pulp, ortoolpy が必要です。環境構築は、組合せ最適化を使おうを参考にしてください。
航空機における搭乗戦略の最適化
論文「航空機における搭乗戦略の最適化」の問題を使わせてもらいましょう。
乗客6人を3グループに分ける。グループ1から順番に搭乗する。同一グループ内は、ランダムな順番とする。総搭乗時間を最小化したい。
考え方
座席$i$の乗客が先に搭乗し、座席$j$の乗客が後に乗った場合の混雑度$a_{ij}$を要素とする行列$A$が与えられるものとします。ここでは、単純に、この混雑度の総和を最小化(=混雑しない組合せの混雑度の最大化)をしましょう。また、同一グループでは、混雑度を半分とします。
混雑しない組合せの混雑度 = (j,iの順番に搭乗する場合の$a_{ij}$の和) + (i,jが同一グループの場合の$a_{ij}$の和)/2
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)かどうかを管理する変数表を作成します。
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]
| Pos | Group | Var | |
|---|---|---|---|
| 0 | 0 | 1 | v1 |
| 1 | 0 | 2 | v2 |
| 2 | 0 | 3 | v3 |
混雑度Aがかかるかどうかを管理する変数表を作成します。
VarNは、座席Firstに後に搭乗し、座席Secondに先に乗った場合を、
VarHは、座席Firstと座席Secondが同一グループの場合を表すものとします。
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]
| First | Second | A | VarN | VarH | |
|---|---|---|---|---|---|
| 0 | 0 | 1 | 0.720 | v19 | v48 |
| 1 | 0 | 3 | 0.302 | v20 | v49 |
| 2 | 0 | 4 | 0.147 | v21 | v50 |
定式化して解きます。制約式の説明は、面倒なので省略。
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] # 解の表示
| Pos | Group | Var | Val | |
|---|---|---|---|---|
| 1 | 0 | 2 | v2 | 1.0 |
| 4 | 1 | 2 | v5 | 1.0 |
| 8 | 2 | 3 | v9 | 1.0 |
| 9 | 3 | 1 | v10 | 1.0 |
| 14 | 4 | 3 | v15 | 1.0 |
| 15 | 5 | 1 | v16 | 1.0 |
座席(Pos)ごとのグループ(Group)が求まりました。
このアプローチは、離散変数が多いので、乗客数が増えると、計算時間が爆発します。その場合は、元の論文のように局所探索法などの近似解法が有効でしょう。
ORセミナーの紹介
11/12(土)に開かれる第3回ORセミナー(Python言語によるビジネスアナリティクス)では、最適化・統計分析・機械学習などのオペレーションズ・リサーチの分野で必要なツールの紹介があります。
このセミナーの参加者特典として、2016年度と2017年度の年会費＋入会費が免除で、学会正会員になれます。正会員になると、上記、機関誌(年12冊)や論文誌を受取れたり、シンポジウムの無料参加などの特典があります。
以上