f:id:ikedaosushi:20190606171352p:plain

PyCon JPの運営メンバーとして自分は、昨年度のPyCon JP 2018のタイムテーブル決めに組み合わせ最適化問題を用いました。ちなみに最終的なタイムテーブルはアルゴリズムで算出された結果を人間がレビューして調整しています。昨年度は時間の制約があり、いくつか反省点があったので今年は改善したいと考えています。

そういうわけで事前調査も兼ねて、カンファレンスのタイムテーブル決めを組み合わせ最適化問題として考え、それをPythonのライブラリを使って解く方法を書きたいと思います。

解きたい問題

カンファレンスのタイムテーブル決めでは下図のように部屋と時間が決まっていて、そこのトークを割り当てていくことが一般的です。

f:id:ikedaosushi:20190602120019p:plain

このとき、部屋の大きさ、時間帯やトークの内容が全く同一であれば何も考えずに隅から配置していけば良いのですが、現実には以下のような条件が発生します。

それぞれのスピーカーが参加できる時間帯(日付)に制限がある
部屋の大きさが異なり、人が集まりそうなトークは大きい部屋に配置したい
同じ時間帯でジャンル被りを少なくしたい
同じ時間帯の発表言語を均等にしたい

こういった条件を加味した上で最適なタイムテーブルを組みたいというのが今回の問題です。カンファレンスのタイムテーブルを例として上げていますが、例えば「ある学校でそれぞれの授業を時間・教室に割り当てる」と言った問題も同じように考えることが可能です。

また、今回はこの条件の1. 2. を加味した問題として解きます。3. 4. を考慮すると線形計画問題でなくなるので、さらに工夫する必要がありボリュームが大きくなりすぎるので今回は線形計画問題で解ける範囲に限定します。

マッチング問題

上記の問題を解くために2部グラフのマッチング問題として扱います。2部グラフのマッチング問題とは今回の例で言うと下記の図のように「トークの集合」と「枠の集合」を「点」で表し「どのトークがどの枠に割り当てられるか」を「辺」で表したグラフを使って「最適な割当」を探す問題のことです。

2部グラフというのは、今回「トーク」と「枠」のように集合が2つ分かれている問題で、これ以降は2部グラフに限定して話を進めます。

f:id:ikedaosushi:20190606152626p:plain

「最適な割当」と書きましたが、「最適」とはなんでしょうか？様々な考え方をすることができるのですが、マッチング問題の中で最も有名な問題として「最大重みマッチング問題」があります。

これは「各辺に重みが与えられているときにその重みを最大化するものを最適としそのようなマッチングを探す」問題です。

例を上げて説明します。次の図を見てください。

f:id:ikedaosushi:20190606155940p:plain

これは、下記のことを表現しています。

「トーク1」は「枠1」には割り当てることができない
「トーク1」は「枠2」「枠3」に割当てることができる
「トーク1」は「枠2」に割当てた方が「重み」が大きくなる

辺の上に書いてある数字が「重み」です。「重み」はそのトークの人気などを考慮して優先度として任意につけることになります。

これを全ての「点」と「辺」に対して考えて、その「重み」の合計が最も大きくなる組み合わせが「最大重みマッチング」になります。

f:id:ikedaosushi:20190606160156p:plain

上記の例のように「重み」の合計が最も大きくなる組み合わせを探す問題が「最大重みマッチング問題」です。

定式化

最大重みマッチング問題の説明をしました。それでは今回解きたい問題を最大重みマッチング問題として定式化して、実際に解いていきます。数式が出てきますが、わからなくても大丈夫なのでそんなに興味ない方は流し読みしてください。

事前準備として次のように定義します。

$G(V, E)$ : トーク(点)の集合と枠(点)の集合を持つ2部グラフ
$w(e)$ : 辺 $e \in E$ の重み
$\delta (v)$ : 点 $v \in V$ に接している辺の集合
$x(e) \in {0,1}$ : 辺 $e \in E$ がマッチングに含まれるとき1をとる変数

このとき下記のように定式化することができます。

$\begin{equation} \begin{aligned} \max_{x} \quad & \sum_{e \in E} w(e)x(e)\\ \textrm{s.t.} \quad & \sum_{e \in \delta(v)}x(e) = 1 \quad \forall v \in V \\ & x(e) \in \{0,1\} \end{aligned} \end{equation}$

この問題を解いていきたいと思います。

Pythonで解く

下記の例を解きます。トーク4つと枠4つ、そしてそれぞれに重みを持っているグラフです。この中で重みが最大になるマッチングを探す問題をGoogle OR-Tools/Pulp/munkresで解きます。

f:id:ikedaosushi:20190606161823p:plain

このグラフを実装上は重み付き隣接行列で表します。上記のグラフは下記の配列になります。通常は0を接続なしとして扱うのですが、今回は0も重みとして扱いたいので接続していない場合はNAを使います。

weights = [
    [6, 7, 15, "NA"],
    [3, 8, 5, 6],
    [12, "NA", 9, 10],
    [4, 11, 9, 11]
]

それでは順番に見ていきましょう。

Google OR-Toolsで解く

Google OR-ToolsはGoogleが公開している最適化問題用のソルバです。PythonやJava、C++などから使うことができます。

developers.google.com

最大重みマッチング問題用のクラスがあるのでそれを使って解くことができます。

from ortools.graph import pywrapgraph

def solve_by_ortools(weights):
    # 最小化問題として解くために最大値から値を引いたものにする
    max_weight = max([max([v for v in row if v != "NA"]) for row in weights])
    weights = [[max_weight - v if v != "NA" else "NA" for v in row ] for row in weights]

    assignment = pywrapgraph.LinearSumAssignment()

    for speaker_idx in range(len(weights)):
        for schedule_idx in range(len(weights[0])):
            weight = weights[speaker_idx][schedule_idx]
            if weight != "NA":
                assignment.AddArcWithCost(speaker_idx, schedule_idx, weight)

    # solve_statusで問題が解けたか判定する
    solve_status = assignment.Solve()
    if solve_status == assignment.OPTIMAL:
        total_weight = len(weights) * max_weight - assignment.OptimalCost()
        assign_weights = {}
        for left in range(0, assignment.NumNodes()):
            right = assignment.RightMate(left)
            assign_weights[(left, right)] = max_weight - assignment.AssignmentCost(left)
        return total_weight, assign_weights
    elif solve_status == assignment.INFEASIBLE:
        raise Exception("No assignment is possible.")
    else:
        raise Exception("Other exception.")

total_weight, assign_weights = solve_by_ortools(weights)
print("合計重み:", total_weight) # => 合計重み: 46
print("割当:", assign_weights) # => 割当: {(0, 2): 15, (1, 1): 8, (2, 0): 12, (3, 3): 11}

割当結果を可視化すると次のようになります。

f:id:ikedaosushi:20190606163658p:plain

Pulp、munkresで解いた結果も同じになるので、後は検算として解いてみたいと思います。

Pulpで解く

Pulpは最もポピュラーなPythonでソルバを扱うライブラリです。Pulp自体がソルバではなく、バックエンドでOSSや商用のソルバを使います。

github.com

Pulpでは次のように解けます。

from pulp import LpProblem, LpMaximize, LpVariable, LpInteger, LpStatus, lpSum, value

def solve_by_pulp(weights):
    from_nodes =  [u for u in range(len(weights))]
    to_nodes =  [v for v in range(len(weights[0]))]

    prob = LpProblem("", LpMaximize)
    choices = LpVariable.dicts("e", (from_nodes, to_nodes), 0, 1, LpInteger)

    objective = []
    for u in from_nodes:
        for v in to_nodes:
            if weights[u][v] != "NA":
                objective.append(weights[u][v] * choices[u][v])

    prob += lpSum(objective), ""

    for u in from_nodes:
        prob += lpSum([choices[u][v] for v in to_nodes]) <= 1

    for v in to_nodes:
        prob += lpSum([choices[u][v] for u in from_nodes]) <= 1

    prob.solve()
    
    solve_status = LpStatus[prob.status]
    if solve_status == "Optimal":
        total_weight = value(prob.objective)
        assign_weights = {}
        for var in prob.variables():
            if var.varValue:
                _, u, v = var.name.split("_")
                u, v = int(u), int(v)
                assign_weights[(u, v)] = weights[u][v]    

        return total_weight, assign_weights
    elif solve_status == "Infeasible":
        raise Exception("No assignment is possible.")
    else:
        raise Exception("Other exception.")

total_weight, assign_weights = solve_by_pulp(weights)
print("合計重み:", total_weight)  # => 合計重み: 46
print("割当:", assign_weights) # => 割当: {(0, 2): 15, (1, 1): 8, (2, 0): 12, (3, 3): 11}

同じ結果が得られました。

munkresで解く

最後はmunkresで解いてみます。おそらくmunkresという名前に耳馴染み無い方も多いと思いますが、ハンガリアン法(こっちは聞いたことがある方もいると思います)というマッチング問題を解くアルゴリズムが実装されたライブラリです。

ハンガリアン法自体はシンプルなアルゴリズムで、ライブラリに頼らなくても実装可能です。自分も学生のときに実装しました。今回はテストされているライブラリを使います。

github.com

en.wikipedia.org

次のように実装することができます。

import sys
from munkres import Munkres, make_cost_matrix, DISALLOWED

def solve_by_munkres(weights):
    # munkresで解くためのデータ形式に変換
    # 同時に最小化問題として解くための処理も行う
    cost_matrix = make_cost_matrix(weights, lambda weight: (sys.maxsize - weight) if weight != "NA" else DISALLOWED)
    m = Munkres()
    indexes = m.compute(cost_matrix)
    total_weight = sum(weights[u][v] for u, v in indexes)
    assign_weights = {(u, v): weights[u][v] for u, v in indexes}
    
    return total_weight, assign_weights

total_weight, assign_weights = solve_by_munkres(weights)
print("合計重み:", total_weight)  # => 合計重み: 46
print("割当:", assign_weights) # => 割当: {(0, 2): 15, (1, 1): 8, (2, 0): 12, (3, 3): 11}

無事、同じ結果を得ることができました。

さいごに

カンファレンスのタイムテーブル決めをマッチング問題として解く考え方と、Google OR-Tools、Pulp、munkresを使って実際にPythonで解く方法を紹介しました。

今年のPyCon JPは最適化問題を使ってタイムテーブル決めを行う予定ですが、もし他のカンファレンスでも需要があれば喜んでご協力するので教えてください。妄想ですが、汎用的なものができたらOSS化したいと思っています。

コードは全てGitHubにあげてあります。

github.com

PyCon JP 2019のチケットはまだ発売されていませんが、送られてきているCFPを見るとかなり盛り上がると思うので、是非チェックしてみてください。運営メンバーも募集中です。組み合わせ最適化の協力も大歓迎です。また昨年に引き続きカンファレンスサイトもNuxt.jsで実装することになりそうなので、そちらも協力してくださる方いたら大歓迎です。

pycon.jp

組み合わせ問題の理論的な部分については下記の書籍を参考にさせていただきました。

amzn.to