Zenn
🤖

Python×Gurobiで遊園地のアトラクションの経路を最適化してみた

2021/08/30に公開
Python
#
networkx
#
matplotlib
#
gurobi
tech

こんにちは。とあるフロントエンドエンジニアです。

学生の頃、研究室のみんなでデ○ズニーランドに行くことになり、「時間内にいっぱいアトラクション乗りたいね〜」「どの順番で当日乗るのが一番効率良いんだろう？？」と話が出ました。

→「これってNP困難やん。じゃあ数理最適化で求めればいっか」と思ったのでサクッと経路を最適化するシステムを作りました。

今回はその定式化部分のプログラムをメモがてらまとめました。

情報系で経路問題に取り組んでいる大学3,4年生の卒業研究の参考にはなるかもしれません笑

同じことをしたいと考えている人も、アルゴリズムは利用できると思うので適宜参考にしてみてください。

プログラム作成手順

システムを実装するにあたり、以下の手順で作成しました。

  1. 技術選定
    • 今回はPythonと数理最適化ソルバーのGurobi Optimizerを使用しました。
    • 使用したアルゴリズムのベースとしては線形計画法の非対称巡回セールスマン問題、Miller-Tucker-Zemlinモデルのポテンシャル定式化パターンを応用。
  2. やりたいことを日本語にまとめる
    • まずはこのシステムで叶えたい内容、制約、どのような経路を最適とするかといった頭の中で描いているものを文字にまとめました。
  3. 定式化
    • 闇雲にプログラムを書いても時間の無駄なのでまずはアルゴリズムを数式に起こしました。
  4. 実装
    • 定式化した数式通りにプログラムを実装
  5. プロット
    • networkxなりmatplotlibなりを使ってネットワークグラフを出力したら完成(これは自己満足のために可視化しただけ)

Gurobi Optimizerとは

膨大な計算量の問題を解いてくれるソルバーです。
他のソルバーとの違いとしては、線形問題の計算に強いとされています。学生であれば大学のライセンスで使えるかもしれません。一度教授に聞いてみましょう。

もう少し詳しい説明はこちら↓
https://www.octobersky.jp/products/gurobi

今回やりたいこと・制約

  1. 開園時間から閉園時間までの時間でより多くのアトラクションに乗りたい。(目的関数: 最大化)
  2. エントランスから出発して最後は必ずエントランスに戻るルートを出したい。
  3. いっぱい乗りたいけど移動距離も考慮したい。→あんまり歩きたくない。
  4. 各アトラクションごとに乗りたい度合いがある。(重み)
  5. 各アトラクションには待ち時間が存在。
  6. アトラクションに乗る = 待ち時間 + サービスタイム
  7. 歩速は一定。

Gurobiのインポート

import gurobipy as gu

以降、本記事でGurobiのメソッドを使用する際には「gu」というモジュールが頻出します。

定式化

やりたいことがまとまったので必要なデータを基に定式化、実装します。

定数

必要なデータをまとめます。
定数の中身は適宜置き換えてください。

# アトラクションの数
NUM_NODE

# 歩く速度(m/m)
# 今回は分速60mで歩く想定
WALK_SPEED = 60

# アトラクション間の距離を格納した辞書
# DISTANCES[0,1] = 100であれば、0から1までは100mを表す。
DISTANCES = {}

# 時間の辞書
# MOVING_TIME[0,1] = 7であれば、0から1まで7分かかることを表す。
MOVING_TIME = {}

# アトラクションを楽しめる時間(サービスタイム)を格納した配列(アトラクションの番号順)
# 言い換えればスプラッシュマウンテンに乗っている大体の時間
SERVICE_TIME = []

# アトラクションの待ち時間を格納した配列(アトラクションの番号順)
WAITING_TIME = []

# 各乗り物の乗りたい度合い(1~5)を格納した配列(アトラクションの番号順)
WEIGHT = []

# 遊園地の開園から閉園までの時間(10:00~21:00)
MAX_PLAY_TIME = 660

# BigM
BIG_M = 10000

モデルの定義

Model = gu.Model("atsp") # モデルの定義

変数

今回は以下の3つの変数を用意します。

x,y,z = {},{},{}
    for i in range(NUM_NODE):
        for j in range(NUM_NODE):
            if i != j:
                x[i,j] = Model.addVar(vtype="B", name=f"枝:x[{i},{j}]") # 枝ijを使用するか否かを表す0-1変数
                y[i,j] = Model.addVar(vtype="C", name=f"時刻：y[{i},{j}]") # アトラクションjの到着時刻
                z[i,j] = Model.addVar(vtype="C", name=f"移動距離:z[{i},{j}]") # アトラクションjまでの移動距離

制約

続いて、先ほど日本語にまとめた制約をプログラムに変えます。

1. 次数制約

まずはネットワークグラフの各枝に関する制約を定式化します。

出発点には1本の枝が入出することを表す制約。
ここからスタートしてここに帰ってくるよ〜という意味。

def degree_constraint(NUM_NODE):
    Model.addConstr(gu.quicksum(x[0,j] for j in range(1,NUM_NODE)) == 1)
    Model.addConstr(gu.quicksum(x[j,0] for j in range(1,NUM_NODE)) == 1)

出発点を除いた各ノードには2つの枝が接続される制約。
各アトラクションには最大で1回訪問できるよ〜(訪問しなくてもいい)という意味。

def degree_constraint_exclude_departure(NUM_NODE):
    for i in range(1,NUM_NODE):
        Model.addConstr(gu.quicksum(x[i,j] - x[j,i] for j in range(NUM_NODE) if j != i) == 0)
        Model.addConstr(gu.quicksum(x[i,j] for j in range(NUM_NODE) if j != i) <= 1)

2. 到着時間に関する制約

ノードjの発時刻の方がノードiの発時刻よりも遅いことを表す制約。
イメージ、アトラクションiを10:00に出発しているのに次に乗るアトラクションjを9:30に出発する解が出てたらそれはヤバイですよね。。。ということを制限する。
次のアトラクションkを出発する時間は、jからkまでの移動時間とkのサービスタイムとkの待ち時間の和で計算します。とりあえず。

def arrival_time_constraint(NUM_NODE, MOVING_TIME, SERVICE_TIME, WAITING_TIME):
    for j in range(1, NUM_NODE):
        Model.addConstr(gu.quicksum(y[i,j] for i in range(NUM_NODE) if i != j) + gu.quicksum( (MOVING_TIME[j,k] + SERVICE_TIME[k] + WAITING_TIME[k]) * x[j,k] - y[j,k] for k in range(NUM_NODE) if k != j) <= 0)

3. 出発点と接続されているノードに関する制約(時間)

出発点の次の点は移動時間+サービスタイム+待ち時間分の時間がかかることを表す制約。
先程との違いは出発点に関連するということ。
もう少し簡単にいうと、出発点での変数yは0にしておきたいという制約。

def constraint_between_departure_2_next_time(NUM_NODE, MOVING_TIME, SERVICE_TIME, WAITING_TIME):
    Model.addConstr(gu.quicksum((MOVING_TIME[0,k] + SERVICE_TIME[k] + WAITING_TIME[k]) * x[0, k] - y[0,k] for k in range(1, NUM_NODE)) <= 0)

4. 移動距離に関する制約

ノードjまでの総移動距離の方が、ノードiまでの総移動距離よりも長いことを表す。
イメージ、アトラクションiまでに合計400m歩いているのに次に乗るアトラクションjまでに歩いた総移動距離が200mの解が出てたらそれはヤバイですよね。。。ということを制限する。

def distance_constraint(NUM_NODE, DISTANCES):
    for j in range(1, NUM_NODE):
        Model.addConstr(gu.quicksum(z[i,j] for i in range(NUM_NODE) if i != j) + gu.quicksum( DISTANCES[j,k] * x[j,k] - z[j,k] for k in range(NUM_NODE) if k != j) <= 0)

5. 出発点と接続されているノードに関する制約(距離)

出発点の次の点には移動距離があることを表す。
こちらは制約3の距離バージョン。
簡単にいうと、出発点での変数zは0にしておきたいという制約。

def constraint_between_departure_2_next_distance(NUM_NODE, DISTANCES):
    Model.addConstr(gu.quicksum(DISTANCES[0,k] * x[0, k] - z[0,k] for k in range(1, NUM_NODE)) <= 0)

6. 非負制約

x[i,j]が1の時だけ、y[i,j],z[i,j]が正の値がとれることを表す制約。
アトラクションiからアトラクションjに行くプランはないのにiからjまで3分歩いた！(変数yに対応)とかいうふざけた解が出ないよう制限する制約。

def nonnegative_constraint(NUM_NODE, BIG_M):
    for i in range(NUM_NODE):
         for j in range(NUM_NODE):
             if i != j:
                 Model.addConstr(z[i,j] <= BIG_M * x[i,j])
                 Model.addConstr(y[i,j] <= BIG_M * x[i,j])

7. 遊ぶ時間の最大を制限する制約

アトラクションの巡回時間の合計は開園から閉園までの時間内に収める。

def play_time_max_limit_constraint(MOVING_TIME, SERVICE_TIME, WAITING_TIME, MAX_PLAY_TIME):
    Model.addConstr(gu.quicksum((MOVING_TIME[i,j]  + SERVICE_TIME[j] + WAITING_TIME[j]) * x[i,j] for (i,j) in x) <= MAX_PLAY_TIME)

目的関数

基本的にはアトラクションの乗りたい度合いを考慮した乗車アトラクション数の最大化。
だがこの場合、最後のアトラクションから閉園時間まで謎の待ち時間が発生してしまう解が出てしまうため、変数yに関するペナルティを課す。
移動距離の変数zのペナルティは、これがないとひたすら歩かされるため、あんまり歩きたくないよ〜という意味。0.5は歩きたくない度合いと思っていただければ。。。

def objective(NUM_NODE, WEIGHT, MOVING_TIME, SERVICE_TIME, WAITING_TIME):
    Model.setObjective(gu.quicksum(WEIGHT[j] * (MOVING_TIME[i,j] + SERVICE_TIME[j] + WAITING_TIME[j]) * x[i,j] for (i,j) in x) - gu.quicksum(y[i,0] for i in range(1,NUM_NODE)) - 0.5 * gu.quicksum(z[i,0] for i in range(1,NUM_NODE)), gu.GRB.MAXIMIZE)

モデルの全体感

def set_model(NUM_NODE, DISTANCES, MOVING_TIME, WALK_SPEED, SERVICE_TIME, WAITING_TIME, WEIGHT):
    Model = gu.Model("tsp")
    x,y,z = {},{},{}
    for i in range(NUM_NODE):
        for j in range(NUM_NODE):
            if i != j:
                x[i,j] = Model.addVar(vtype="B", name=f"枝:x[{i},{j}]") # 枝ijを使用するか否かを表す0-1変数
                y[i,j] = Model.addVar(vtype="C", name=f"時刻：y[{i},{j}]") # アトラクションjの到着時刻
                z[i,j] = Model.addVar(vtype="C", name=f"移動距離:z[{i},{j}]") # アトラクションjまでの移動距離
    Model.update()

    # 1 次数制約
    degree_constraint(NUM_NODE)
    degree_constraint_exclude_departure(NUM_NODE)
    # 2 到着時間に関する制約
    arrival_time_constraint(NUM_NODE, MOVING_TIME, SERVICE_TIME, WAITING_TIME)
    # 3 出発点と接続されているノードに関する制約(時間)
    constraint_between_departure_2_next_time(NUM_NODE, MOVING_TIME, SERVICE_TIME, WAITING_TIME)
    # 4 移動距離に関する制約
    distance_constraint(NUM_NODE, DISTANCES)
    # 5 出発点と接続されているノードに関する制約(距離)
    constraint_between_departure_2_next_distance(NUM_NODE, DISTANCES)
    # 6 非負制約
    nonnegative_constraint(NUM_NODE, BIG_M)
    # 7 遊ぶ時間の最大を制限する制約
    play_time_max_limit_constraint(MOVING_TIME, SERVICE_TIME, WAITING_TIME, MAX_PLAY_TIME)

    #目的関数
    objective(NUM_NODE, WEIGHT, MOVING_TIME, SERVICE_TIME, WAITING_TIME)

    Model.update()
    Model.__data = x,y,z
    return Model

この関数で定数、変数、制約を定義したモデルをreturnします。

最適化

最後は最適化するためのコードを書いて終了です。

if __name__ == "__main__":
    Model = set_model(NUM_NODE, DISTANCES, MOVING_TIME, WALK_SPEED, SERVICE_TIME, WAITING_TIME, WEIGHT)
    Model.optimize()

求められた解をmatplotlibやらnetworkxやらで可視化してあげるとこんな感じでプロットできます。(例)↓
最適化のネットワークグラフ例

著作権の都合上載せられませんが、当時私は遊園地のマップと上のような出力画像をレイヤードして一緒に出力するようにしました。

アトラクションなどの座標の取得は緯度経度でも良いのですが、園内マップのエントランスからの位置ベクトルで計算するとわりと楽でした。

応用可能なシチュエーション

  • デ○ズニーシー
  • 富○急ハイランド
  • U○Jなど

大抵の遊園地で応用できると思います。
適宜シチュエーションに合った制約を追加して使ってください。

まとめ

正直、時間と距離は比例する部分もあるので遊ぶ時間の最大化だけで良いのでは？と思う部分もありますがそこは今は気にしない。。。笑

以上です。

似たような定式化を考えている際は参考にしていただけると。。。

Discussion