第6回：OR-Toolsで巡回セールスマン問題を解く ～京都弾丸観光ツアーの作成を事例に～【ブレインパッドの数理最適化ブログ】

はじめに

本記事では、京都の観光地を巡る巡回セールスマン問題を例としながら、OR-Toolsという最適化ツールをご紹介します。

OR-Toolsは Google が開発しているフリーの最適化ツールで、幅広い最適化問題を扱えます。中でも今回はメタヒューリスティクスを用いてルーティング問題を解く、routingモジュールを使用します。

https://developers.google.com/optimization/routing

最適化のツールは色々ありますが、OR-Toolsの利点として以下の3つがあげられます。

• 1つ目に、オープンソースのソフトウェアで、無料で使えます。
• 2つ目に、実務での大規模問題にも耐えうる十分な性能があります。最適化のソルバーは有料のもののほうが性能が良いことが一般的ですが、昨年行われたKaggleの最適化コンペの2位の解法に使われたり、solverのコンペで上位に入賞したりといった実績があります。
• 3つ目に、多様な制約を柔軟に入れられます。これは（routing問題を解く際に）裏でconstraint programmingの手法が用いられていることにもよります。

一方でメタヒューリスティクスを用いているため、必ずしも厳密解を求めているわけではないことには注意が必要です（メタヒューリスティクスについては大阪大学の梅谷先生のこちらのスライド（大規模な組合せ最適化問題に対する発見的解法 | PPT）が参考になります）。

なかなか便利なツールだと思いますが、日本語の情報はあまりないので、その点についてもこの記事が参考になればと思います。

さて、本記事では例として、観光地を巡る順番を求める問題を解いていきます。シンプルな巡回セールスマン問題（TSP）は学問として数十年議論されており、私のような文系人間には口を出す余地はありません*1 。この記事では、少しひねりを加えた制約を入れたTSPを解いていきます。

我々が解くようなビジネス上の問題は、複雑な業務上の制約を満たしながら解を求めることが要求されます。ビジネスに数理最適化を取り入れる困難さは、①そもそも何をどう最適化問題に落とし込むか、②複雑な制約をどう数式やコードで表現するか、の2つにあると思っています。前者についてはブレインパッドの数理最適化連載の第4回や梅谷先生の社内講演会レポートをぜひご覧いただければと思います。本記事で具体例を解くことを通じて、後者について、つまり弊社の分析官が「現実の問題で出てくる制約に対して、どう工夫して対処しているか」についての紹介にもなれば、というのが裏テーマです。

問題設定

具体的な問題として、京都の観光地を1日でできるだけたくさん巡る問題を解いていきます。どの観光地に行くか、それらをどのような順番で回るか、を決める問題です*2。

以下のような制約を置いて解いていきたいと思います。徒歩でガンガン観光地を回る弾丸ツアーになります笑

• なるべくたくさんの観光地を見ることを目的とする
• スタートは朝5時半に京都駅出発（たぶん東京発の夜行バスが京都駅に着くのがそのくらいの時間）
• ゴールは夜8時に京都駅到着
• 移動手段は徒歩のみ。歩くのにも楽しい街ですからね（夏と冬を除く）
• 移動時間だけでなく、観光する時間も考慮する
• 観光地が開いている時間に行かなければならない

以上の条件で解いてみた後、さらにいくつか制約を追加して解いていきます。

データの取得・準備

京都市のオープンデータより、観光地のデータをとってきます*3。このデータには、各観光地について

• 名称
• 住所
• 開始/終了時刻
• 想定所要時間

などの情報があります（正確に言うと、csvに列名の記載を見つけられなかったので、「おそらくこの列はこういう意味であろう」という勘で使用しています。したがってこの情報が正しい保証はありません。また、更新がしばらくされておらず若干古いデータもあるようなので、実際に行く際は最新の情報をご確認ください。）

位置情報（緯度・経度）はもとのデータにはなかったので、Google Maps APIを用いて付与しました。情報が一部欠けていた観光地は除いています。また、開始地点として、京都駅を手で追加しています（結果37か所になりました）。

整形したものが以下のようなデータになります。なお、時間は時刻表示でなく分単位に変換しておくと最適化の計算の際に便利です。

地図に描くとこんな感じです。

最適化に用いる際に使いやすいように、辞書として必要な情報を取り出しておきます

required_time = places.set_index('name')['required_time'].to_dict()
open_time_min = places.set_index('name')['open_time_min'].to_dict()
close_time_min = places.set_index('name')['close_time_min'].to_dict()
list_name = places['name'].to_list()

また、任意の2地点間の移動時間も作成して辞書に格納します。京都はだいたい碁盤状に道があるので、manhattan距離を用いて計算しています。

def calc_manhattan_dist(lat1, lon1, lat2, lon2):
    # 北緯35度での経度1度のメートル91287.7885、緯度1度のメートル110940.5844を用いる
    dist_meter_x = np.abs(lon1 - lon2) * 91287.7885
    dist_meter_y = np.abs(lat1 - lat2) * 110940.5844
    return dist_meter_x + dist_meter_y

def calc_time_matrix(places):
    walk_time_matrix = dict()
    for f_index, f_row in places.iterrows():
        for t_index, t_row in places.iterrows():
            dist_meter = calc_manhattan_dist(f_row.lat, f_row.lon, t_row.lat, t_row.lon)
            walk_time_matrix[(f_row['name'], t_row['name'])] = (np.ceil(15 * dist_meter / 1000))
    return walk_time_matrix

walk_time_matrix = calc_time_matrix(places)

OR-Toolsによる最適化

それでは、実際に上の問題を解いていきます。以下のステップで説明していきます。

1. index managerの作成
2. モデルの記述
3. 求解
4. 解の確認
5. 制約を追加して再び求解

ということで、一旦上記問題設定で解いた上で、さらにいろいろな制約を追加して解いてみます。

なお、OR-Toolsについては

• 公式のガイド（https://developers.google.com/optimization/routing）
• API documentation（http://google.github.io/or-tools/python/ortools/constraint_solver/pywrapcp.html）

もご参照いただければと思います。

Step1：index manager/node managerの作成

OR-Toolsのrouting modelでは、巡るべき地点を表すものにNode/Indexのふたつがあります。Nodeは我々が認識している地点番号、Indexはソルバーの内部で処理される際に用いられるものです。ややこしいですが、ソルバーの内部で地点が複製されたりもするため、このように使い分ける必要があります。NodeとIndexの紐づけを管理してくれるのがIndexManagerです。また、（移動時間などの）定数を取得する際には地点番号よりも地点名（Name）のほうが扱いやすいため、IndexManagerと同様の働きを持つNodeManagerを作成しています。

RoutingIndexManagerを作成する際には地点数、vehicle数、start/end地点のNodeを指定します。

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

class NodeManager:
    """Node(int)とName(文字列)を紐づける処理を担う"""
    def __init__(self, list_name):
        self.list_name = list_name
        self.n = len(self.list_name)
        self.list_node = list(range(self.n))

        self.name_map_node = dict(zip(self.list_name, self.list_node))
        self.node_map_name = {v: k for k, v in self.name_map_node.items()}

    def NodeToName(self, node):
        return self.node_map_name[node]

    def NameToNode(self, name):
        return self.name_map_node[name]

# node managerの作成
node_manager = NodeManager(list_name)
# index managerの作成
index_manager = pywrapcp.RoutingIndexManager(
    node_manager.n, 1, node_manager.NameToNode('京都駅'))

Step2：モデルの記述

IndexManagerインスタンスを引数にしてroutingモデルを作成します。ここに問題の色々な設定を記述していきます。

# routingモデルの作成
routing = pywrapcp.RoutingModel(index_manager)

# 移動時間callbackの作成
def time_callback(from_index, to_index):
    # ソルバーからもらったIndexを、移動時間行列にアクセスできるようNameに変換する
    from_node = index_manager.IndexToNode(from_index)
    to_node = index_manager.IndexToNode(to_index)
    from_name = node_manager.NodeToName(from_node)
    to_name = node_manager.NodeToName(to_node)

    # 移動時間 + 観光地の観光時間が次の地点の到着までに加算される
    return walk_time_matrix[from_name, to_name] + required_time[from_name]

transit_callback_index = routing.RegisterTransitCallback(time_callback)
# Arcコストの設定
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

# 観光地に行かない場合はペナルティがかかる
for i in node_manager.list_name:
    if i == '京都駅':
        continue
    index = index_manager.NodeToIndex(node_manager.NameToNode(i))
    routing.AddDisjunction([index], 10000)

# 時間Dimensionを作成
routing.AddDimension(
    transit_callback_index,
    60,  # 待ち時間を最大60分まで許容する
    60 * 20,  # 最大時間（=Endに帰って来る期限）
    False,  # 時間は0から始まる必要はない
    'Time')  # Dimensionの名前
time_dimension = routing.GetDimensionOrDie('Time')

# 営業時間枠に収まる制約
for i in node_manager.list_name:
    index = index_manager.NodeToIndex(node_manager.NameToNode(i))
    open_time = open_time_min[i]
    close_time = close_time_min[i]
    time_dimension.CumulVar(index).SetRange(int(open_time), int(close_time))

地点間の移動にかかるコストを設定するのに、(1) indexを引数、変化量を出力とする関数を作成、(2) callbackとして登録、(3) ArcCostに設定、する必要があります。ArcはグラフのEdgeに対応します。すなわち、ここでは移動時間の総和をコスト（最適化の目的関数）として設定しています。なお、移動時間callbackを作成する際に、徒歩での移動時間に加えてその観光地の観光にかかる所要時間も考慮します。

今回はなるべく多くの観光地を巡ることを目的としますが、そのためには近い観光地を効率よく回る必要があるため、近さの概念を目的関数に組み込むことで効率的に解くことができます。

AddDisjunction

今回の問題では、すべての観光地に行くことはできません。沢山ある観光地の候補の中から、なるべく多く回れるように行く地点を選択する必要があります。AddDisjunctionを用いることで、「その地点には行かなくてもいい」と設定できます*4。

このとき最適化の目的関数は、「先程設定した移動時間の総和+行かなかった地点数×ペナルティ」となります。ペナルティが小さすぎると、どこも訪れないのが最適解になってしまうので、十分大きな数（ここでは10000）をペナルティに設定しておきます。

Dimension

次に、時間枠内に各地点に訪れる制約を課します。OR-Toolsでは、時間や積載量などの地点を経るごとに蓄積される変数を追うのにDimensionを用います。先ほど登録したcallbackを用いてDimensionを作成します。

Dimensionの各地点での値に制約を入れたいときは、dimension_name.CumulVar(index)でアクセスできます。例えば時間枠に収まる制約を入れるには、SetRangeを用いて、ここからここの間に到着時間が入る、という制約を表現できます。なお、RoutingModelは整数しか受け付けないものがあるためint型にしていることにご注意ください。

Step3：求解

問題の記述が終わったら、serach_parameterを決めた上で求解します。

# 解法の設定
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.time_limit.seconds = 10
search_parameters.local_search_metaheuristic = (
    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)

# 最適化の実施
solution = routing.SolveWithParameters(search_parameters)

serach_parameterは基本的にデフォルトのままでよいと思いますが、Disjunctionが用いられていると挙動が変になることがあるという話もある（ https://activimetrics.com/blog/ortools/exploring_disjunctions/）ので、一応メタヒューリスティクスとしてGUIDED_LOCAL_SEARCHを選択しています。GUIDED_LOCAL_SEARCHを用いると、（たとえ良い解にたどり着いていたとしても）解の探索がなかなか終わらないため、時間の上限を入れるのがベターです。

なお、使用できるメタヒューリスティクスの一覧はこちら、用いられているデフォルトのヒューリスティクスは（やや見ずらいですが）こちらで確認できます。

Step4：解の確認

def print_solution(routing, solution):
    """結果を出力"""
    vehicle_id = 0
    index = routing.Start(vehicle_id)
    while not routing.IsEnd(index):
        previous_index = index
        name = node_manager.NodeToName(index_manager.IndexToNode(previous_index))
        output = name + ': '
        output += minutes2hm(solution.Value(time_dimension.CumulVar(index)))
        print(output)
        index = solution.Value(routing.NextVar(index))
    print('京都駅:', minutes2hm(solution.Value(time_dimension.CumulVar(index))))
    
print_solution(routing, solution)

# 出力は以下
# 京都駅: 05:30
# 下鴨神社: 07:09
# 北野天満宮: 08:37
# 地蔵院: 09:22
# 雨宝院: 10:05
# 妙蓮寺: 10:40
# 宝鏡寺: 11:16
# 相国寺: 12:01
# 廬山寺: 13:24
# 安養寺: 14:27
# 禅居庵: 14:56
# 南岩倉明王院不動寺: 15:17
# 蓮光寺: 15:31
# 矢田寺: 16:29
# 清水寺: 17:54
# 京都駅: 19:24

solution.Value(x)で、解におけるxの値を取得できます。NextVar(Index)がindexの次に行く場所を示すので、while文でEnd地点に行き着くまでループすることで、解における一連のルートを取得できます。Dimensionの値は、制約を入れたときと同じようにdimension_name.CumulVar(index)でアクセスします。

なお、routing.Start(vehicle_index)は解に依存しない（どんな解でもstartは京都駅である）ため、solution.Valueにつっこむ必要はありません。

地図上にプロットすると以下のような結果になります。下の方にある青い丸が京都駅です。ぱっと見では少し違和感があるルートかもしれませんが、朝イチだと開いてない場所も多いため、朝早くから開いている下鴨神社から攻めています。

Step5：制約の追加

以上で基本的な最適化問題を解くことができました。ここからが本題です。より実践的な制約の例として、更に以下の制約を加えて解いていきます。

• 任意の地点で終了
• 昼食をとる
• 世界遺産は朝イチに行かない

ここらへんからストーリーとしては若干無理めな設定が入ってきますが、そこはスルーして続けます笑

任意の地点で終了

今まではスタートもエンドも京都駅、つまり日帰り（あるいは京都駅付近で宿泊）という設定にしていました。これを、エンド地点はどこでもいいものと設定します。想定としては、翌日は古い友人に会う予定があるからどこかの安宿で１泊する、くらいでしょうか。この設定では、京都駅に時間までに帰ることを想定してルートを組む必要はなくて、一方通行で最適なルートを探すことになります。

こちらの制約は、ゴール地点に帰る際の距離を0という風にtime_callbackを書き換えれば対応できます。

def time_callback(from_index, to_index):
    from_node = index_manager.IndexToNode(from_index)
    to_node = index_manager.IndexToNode(to_index)
    from_name = node_manager.NodeToName(from_node)
    to_name = node_manager.NodeToName(to_node)

    if to_name == '京都駅':
        return 0
    return walk_time_matrix[from_name, to_name] + required_time[from_name]

昼食をとる

予めサーチしておいた良さげなお店が、相国寺、安養寺、清水寺の近くにあるものとしましょう。この中のいずれかで昼食の時間を60分取りたいとします*5。

該当の3地点のどこかで追加で60分過ごす、という制約を入れることになりますが、どこかひとつの地点で時間を追加、という制約をどのように表現したらいいでしょうか。

このような制約を入れるときの常套手段に、地点を複製する。というものがあります。3地点について、地点を複製したものをplacesに追加します。名前を変更し、所要時間を1時間、時間枠を11〜14時くらいにしておきましょう。複製した3地点について、「どれかひとつにだけ行く」という制約を加えることで制約を表現できます。

# 昼食候補地点を複製
list_lunch_candidate = ['相国寺', '安養寺', '清水寺', ]
places_lunch = places[places['name'].isin(list_lunch_candidate)].copy()
places_lunch['name'] = places_lunch['name'].apply(lambda x: x + '_lunch')
places_lunch['required_time'] = 60
places_lunch['open_time_min'] = hm2minutes('11:00')
places_lunch['close_time_min'] = hm2minutes('14:00')
places = pd.concat([places, places_lunch], axis=0)

この「どれかひとつにだけ行く」という制約を入れるのに、ActiveVarを使います。ActiveVarはそのindexが訪れられる場合は1を、訪れられない場合は0をとります。3地点のActiveVarの合計がジャスト1になる、とすれば制約を表現できます。

ちなみに、このような数式で表せる制約はAddConstraintで直接入れることができます。

# 昼食候補地点のindexを取得
list_lunch_index = [index_manager.NodeToIndex(node_manager.NameToNode(i)) 
                    for i in list_lunch_candidate]

# 3地点のうちひとつのみ訪れる制約を入れる
condition = (routing.ActiveVar(list_lunch_index[0]) + 
             routing.ActiveVar(list_lunch_index[1]) + 
             routing.ActiveVar(list_lunch_index[2]) == 1)
routing.solver().AddConstraint(condition)

世界遺産は朝イチには行きたくない

これも謎設定ですが、朝イチだとまだ目が覚めていないのですこし時間が経ってから行きたい、ということにしましょう。2通りの方法を紹介します。

list_heritage = ['北野天満宮', '下鴨神社', '仁和寺', '清水寺', ]

# 方法１：counting dimensionを使用    
def count_callback(index):
    return 1

# counting dimensionの作成
count_callback_index = routing.RegisterUnaryTransitCallback(count_callback)
routing.AddDimensionWithVehicleCapacity(
    count_callback_index,
    0,  # null capacity slack
    [100],  # maximum count
    True,  # start cumul to zero
    'count')
count_dimension = routing.GetDimensionOrDie('count')

# 該当する各地点が2番目以降に訪れる制約を入れる
for h in list_heritage:
    h_index = index_manager.NodeToIndex(node_manager.NameToNode(h))
    cond = count_dimension.CumulVar(h_index) >= 2
    routing.solver().AddConstraint(cond)

# 方法２：RemoveValueを使用
start_index = index_manager.NodeToIndex(node_manager.NameToNode('京都駅'))
for h in list_heritage:
    h_index = index_manager.NodeToIndex(node_manager.NameToNode(h))
    routing.NextVar(start_index).RemoveValue(h_index)

ひとつは、counting dimensionを使う方法です。つまり、何番目に着いたかを表すDimensionを作成し、それに対して「世界遺産の地点はcounting dimensionが2以上」という制約を入れます。Dimension作成の際に、time dimensionと違ってfrom/to双方の地点に依存することはないため、AddDimensionWithVehicleCapacityを使用します。

ちなみにこのcounting dimensionも汎用性が高くて、例えば「地点Aと地点Bは3つ以上離れてほしくない」といった制約もかんたんに入れることができます。こちらのブログで色々紹介されていて非常に参考になります。

もうひとつは、NextVarとRemoveValueを使って直接、ここからここへは行くかないように設定する方法です。今回は朝一番で行きたくないということなので、京都駅から該当の観光地に行くことを避ければよいです。

以上3つの制約を入れて解いた結果がこちらになります。

制約に従って

• 最初に行く地点が変わった
• 最後も駅から遠くにある

ことがわかります。また、（この図からは見にくいですが）お昼はちゃんと行ってることも確認できました。結果は14箇所だったので、1日あれば徒歩でもかなり回れるんだなあという印象ですね。

むすび

本記事ではOR-Toollsでルーティング問題を解く方法、およびさまざまな制約を入れる方法について、具体例をもとに紹介させていただきました。OR-Tools特有の方法もありましたが、ノードの複製や、移動時間を用いて制御する方法などは汎用性が高い考え方だと思います。

実際の案件では求解の時間なども含めて、さらにいろいろな制約があります。冒頭に「複雑な制約を入れることが実務に最適化を応用することの困難のひとつある」と書きましたが、これは工夫のしがいがあって非常に面白いポイントでもあります。ブレインパッドでは、共にこのような問題に取り組む仲間を募集しています。興味のあるエンジニア・データサイエンティストのみなさま、ぜひご応募をお待ちしています！！

https://www.brainpad.co.jp/recruit/

*1:ちなみにTSPは解くのが非常に難しい問題とされますが、最近の研究では億単位の地点数の問題が解かれているそうです（Star Tours）。ごいすー

*2:余談ですが、私は大学時代京都に住んでいたわりに、住んでいるといつでも行けると思って、名所にすらろくに行ったことがありませんでした。せっかく社会人になって余裕ができたので、久々に行きたいと思っているのですが、コロナがあれなので、まだしばらく先になりそうです。ただ旅行ってのはプランを妄想するだけでも楽しいですからね。というのが、この問題を解こうと思ったきっかけになります。

と、熱い想いを語ったところですが、ここで残念なお知らせです。実はすでにこの手の問題はやられており、公式の観光サイトがその機能を持っています（ https://ja.kyoto.travel/jtp_map/）。ぐぬぬ。だがしかし、公式のものはシンプルに順番を教えてくれるだけで、そのルートの実現可能性などは考慮されていません。より実践的なプランを妄想する際には色々な制約を加える必要があるので、独自に作る意義はあるでしょう。

*3:出典：「京都観光見もの情報」（京都市 産業観光局 観光ＭＩＣＥ推進室） （https://data.city.kyoto.lg.jp/node/92390）。

*4:RoutingModelは「制約を満たし、かつすべての地点を巡るようなルート」が見つけられない場合は、Noneを出力します。エラーも出ないので注意が必要です。したがって今回のようなケースはAddDisjunctionが必須です

*5:相国寺は有名大学の近く、安養寺は繁華街のあたり、清水寺は言わずもがななので、おそらく近くには美味しいものを食べられる場所が沢山あるだろうということで3つが選ばれました。