Fusic Tech Blog

Fusion of Society, IT and Culture

時刻によって移動時間が変わる巡回セールスマン問題をFixstars Amplify を使って解く
2021/10/25

時刻によって移動時間が変わる巡回セールスマン問題をFixstars Amplify を使って解く

TSP (Travelling Salesmen Problem)は、セールスマンが複数地点を1回だけ通る時に移動距離が最小となる経路を求める問題です。考える変数を距離から時間に変えると、複数地点を通る際の移動時間を最小化することもできます。しかし、この問題は、移動時間が時間によって変わることは想定されていません。今回は、TSP + 時間によって変動する移動時間 という問題を設定し、解いてみます。

イメージしているのはゴミ収集車。ゴミを回収するためにいろいろな地点を回る必要がありますが、渋滞に巻き込まれると大変なので、渋滞情報を加味して経路を決める必要があります。

似たような問題は、VRP(Vehicle Routing Problem) として知られています。VRP は最適化問題の一つで、複数の車両で複数の地点を回る際の移動距離を最小化させる問題です。(参考: https://qiita.com/r_nsd/items/19dcb30f5478384f90d3)

なお、時間情報を加味するというアイディアは、https://arxiv.org/abs/1903.06322 の論文を参考にしました。この論文自体は、CVRP(容量制約付きの運搬経路問題) を対象に書かれていますが、今回は時間の取り扱い方を参照するにとどまっています。将来的には、VRP・CVRP も試してみたいです。

また、時間情報を入れるというところを重要視したので、現状車は1台だけを想定しています。

準備

準備として、各種import とToken の取得を行います。 実行は、Google Colaboratory で行いました。

!pip install amplify
from amplify import BinaryPoly, BinaryQuadraticModel, gen_symbols, sum_poly, Solver, FixstarsClient, decode_solution
from amplify.constraint import equal_to, less_equal, clamp, penalty
import json
import numpy as np

# Token の取得
with open("drive/MyDrive/quantum_tokens/fixstars.txt", "r") as f:
    token = f.read()

token は、https://amplify.fixstars.com/ja/ にサインアップして、取得しました。

前提条件

時間 の概念をイジングモデルに組み込むために、時間を時刻という概念に変換します。

例えば、時刻の刻みの最小値が5分で、A 地点からB 地点まで3分かかる場合、A地点からB地点は5分かかるという取り扱いをします。ここからは、移動時間は、何時刻かかるか? に着目して計算を行います。

3地点に関するVRP を解く場合、以下のようなテーブルを時刻分準備します。

この表の読み方は、T=0 の場合、地点A から地点B に行くのに1時刻、地点Bから地点A に行くのに2時刻 という風に読みます。列方向が到着地点、行方向が出発地点を表しています。

T = 0

A B C
A 0 1 2
B 2 0 2
C 3 3 0

T = 1

A B C
A 0 3 3
B 2 0 2
C 1 2 0

・・・

今回は、とりあえずこんな感じで定義しました。

# かかる時間の定義
ary = np.array([
       [
              [0, 1, 2],
              [2, 0, 2],
              [3, 3, 0],
       ],
       [
              [0, 3, 3],
              [2, 0, 2],
              [1, 2, 0],
       ],
       [
              [0, 3, 3],
              [3, 0, 3],
              [1, 2, 0],
       ],
       [
              [0, 3, 3],
              [3, 0, 3],
              [1, 1, 0],
       ],
       [
              [0, 3, 3],
              [3, 0, 3],
              [1, 1, 0],
       ],
       [
              [0, 3, 3],
              [3, 0, 3],
              [1, 1, 0],
       ],
       [
              [5, 5, 5],
              [5, 5, 5],
              [5, 5, 5],
       ],
       [
              [5, 5, 5],
              [5, 5, 5],
              [5, 5, 5],
       ],
       [
              [5, 5, 5],
              [5, 5, 5],
              [5, 5, 5],
       ],
       [
              [5, 5, 5],
              [5, 5, 5],
              [5, 5, 5],
       ],
       [
              [5, 5, 5],
              [5, 5, 5],
              [5, 5, 5],
       ],
])

最後の方で5時刻を連発しているのは、後述の、T の次元を確保するための作戦です。

変数の定義

それでは、変数を定義していきます。

変数V は、 T x N x N で表されます。

T: 時刻数
N: 巡る必要がある地点の数
V = T x N x N

変数は、

V(T=0)

A B C
Step1 q0 q1 q2
Step2 q3 q4 q5
Step3 q6 q7 q8

V(T=1)

A B C
Step1 q9 q10 q11
Step2 q12 q13 q14
Step3 q15 q16 q17

という感じになっており、何番目にどこに行くか?を表しています。

実際に変数を定義するコードは、以下のとおりです。

N = len(ary[0])
T = len(ary)

variables_n = gen_symbols(BinaryPoly, T, N, N)

制約条件

制約条件はいっぱいあります。

同一時刻に2拠点以上は存在できない

less_than_1_for_a_time = sum([less_equal(sum_poly(t), 1, f'less_than_1_for_a_time===>{{"time": {t}}}===>') for t in variables_n])

1拠点への訪問は1回限り

visit_a_single_place_only_once = sum([
    equal_to(variables_n[:, :, p], 1, f'single_place_only_once===>{{"place": {p}}}===>') for p in range(N)
])

Step1 ~ Step3 は、全て埋まる必要がある

need_to_fill_nth_place = sum([
    equal_to(variables_n[:, p, :], 1, f'nth_place===>{{"place": {p}}}===>') for p in range(N)
])

Step1 ~ Step3は、順番にうまる必要がある

# Step の順番を守る制約条件が必要
step_order = BinaryConstraints()
for step in range(N-1):
    for t in range(T):
        this_variable = sum_poly(variables_n[t, step])
        tmp_current = variables_n[:t+1, step]
        tmp_future = variables_n[:t+1, step+1]
        step_order += clamp(sum_poly(tmp_current) - sum_poly(tmp_future), 0, 1, f'step_order===>{{"step": {step}, "time": {t}}}===>')

移動に必要な時刻を確保する

ここが一番やっかい、かつ、面白いところです。

時刻0, Step1, 地点A が1 であり、地点A -> 地点C にいくのに 3 時刻かかる場合、 3時刻の間 variables_n[t, 1:, 2] は、0にならなくてはいけません。この条件を、penalty 関数を使って表現しています。

avoid_time_overrap = []
for t in range(T):
    for step in range(N):
        if step + 1 == N:
            continue
        for location in range(N):
            location_at_step = variables_n[t, step, location]
            for from_location in range(N):
                # 自分からの移動を考えるので。
                if location != from_location:
                    continue
                for to_location in range(N):
                    # 同じ地点の移動はないので。
                    if to_location == from_location:
                        continue
                    # 時刻t にfrom_location から to_location に行くのにかかる時間
                    time_to_arrive = ary[t][from_location][to_location]
                    for time_consumed in range(time_to_arrive):
                        if t + time_consumed >= T:
                            continue
                        avoid_time_overrap += penalty(
                            location_at_step * variables_n[t + time_consumed, step + 1, to_location], 
                            f'overrap_penalty===>{{"t": {t}, "step": {step}, "location": {location}, "from_location": {from_location}, "to_location": {to_location}, "time_to_arrive": {time_to_arrive} }}===>'
                        )

コスト関数

コスト関数は、最後のStep に到着する時刻がができるだけ早くなるように設定します。

cost = sum([sum_poly(variables_n[t, N-1, :]) * t for t in range(T)])

これで、準備が整いました。

SOLVE!!!

constraints = BinaryConstraints()
constraints += avoid_time_overrap
constraints += step_order
constraints += less_than_1_for_a_time
constraints += visit_a_single_place_only_once
constraints += need_to_fill_nth_place
equation = constraints + 0.001 * cost


equation = BinaryQuadraticModel(sum(avoid_time_overrap) + sum(less_than_1_for_a_time) + sum(visit_a_single_place_only_once) + sum(need_to_fill_nth_place) + sum(step_order) + 0.1 * cost)

timeout = 1

client = FixstarsClient()
client.token = token
client.parameters.timeout = timeout * 1000  # タイムアウト1秒

solver = Solver(client)
result = solver.solve(equation)
if len(result.solution) >= 1:
    print(decode_solution(variables_n, result.solutions[0].values))

問題を解くと、

array([[[1., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 1., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 1.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]],

       [[0., 0., 0.],
        [0., 0., 0.],
        [0., 0., 0.]]])

のような結果が返ってきました。 結果を見やすく整理すると、

answer.sum(axis=0)
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

となります。

この結果から、1番目に A地点、2番目にB地点、3番目にC 地点を訪れれば良いことがわかります。

結果の検証

今回の問題はとてもシンプルなので、全探索で結果を検証してみたいと思います。 少しコードが雑ですが、ご容赦いただけると。。。

考えられる6通りのすべての組み合わせ (以下)で、どれだけ到着までの時間がかかるかを計算しています。

  • ABC
  • ACB
  • BAC
  • BCA
  • CAB
  • CBA
for t in range(T):
    ## A
    # A -> B or A -> C
    ab = ary[t][0][1]
    if ab >= T:
        continue
    print("T = ", t,"abc = ", ary[ab][1][2] + ab + t)
    ac = ary[t][0][2]
    if ac >= T:
        continue
    print("T = ", t,'acb = ', ary[ac][2][1] + ac + t)

    ## B
    # B -> A or B -> C
    ba = ary[t][1][0]
    if ba >= T:
        continue
    print("T = ", t,"bac = ", ary[ba][0][2] + ba + t)
    bc = ary[t][1][2]
    if bc >= T:
        continue
    print("T = ", t,'bca = ', ary[bc][2][0] + bc + t)

    ## C
    # C -> A or C -> B
    ca = ary[t][2][0]
    if ca >= T:
        continue
    print("T = ", t,"cab = ", ary[ca][0][1] + ca + t)
    cb = ary[t][1][2]
    if cb >= T:
        continue
    print("T = ", t,'cba = ', ary[cb][1][0] + cb + t)

結果は、

T =  0 abc =  3
T =  0 acb =  4
T =  0 bac =  5
T =  0 bca =  3
T =  0 cab =  6
T =  0 cba =  5
T =  1 abc =  7
T =  1 acb =  5
T =  1 bac =  6
T =  1 bca =  4
T =  1 cab =  5
T =  1 cba =  6
・・・

という形になりました。

これより、 Aスタートのときは、 A -> B -> C、Bスタートのときは、 B -> C -> A、Cスタートのときは C->B->A の順に回れば良いことがわかります。

Amplify で計算した結果も、 A -> B -> C となっており、正解していることがわかります。

出発地点固定

出発地点を固定するためには、変数を定義した variables_n = gen_symbols(BinaryPoly, T, N, N) のあとに、

variables[n, 0, 1] = BinaryPoly(1)

という形で、1にしたい値を直接代入してあげれば良いです。

この場合の結果は、

answer = decode_solution(variables_n, result.solutions[0].values)
answer.sum(axis=0)
# => array([[0., 1., 0.],
#      [0., 0., 1.],
#      [1., 0., 0.]])

となり、こちらも、上記の答えとぴったりあいます。

まとめ

今回は、少しむずかしいTSP・少し簡単なVRPを、Fixstars Amplify を使って解いてみました。

どの程度大規模な問題が解けるか?はやってみないとわからないところではありますが、解けることはわかってよかったです。

今後は、

  • 問題の大規模化
  • 車複数台対応
  • CVRP 対応

を試していきたいと思います。

あと、ダミーデータを準備するのが大変なことがわかりました。

ちなみに、

from math import factorial
for i in range(3, 21):
    print(i, factorial(i))
3 6
4 24
5 120
6 720
7 5040
8 40320
9 362880
10 3628800
11 39916800
12 479001600
13 6227020800
14 87178291200
15 1307674368000

ということで、15地点の組み合わせは1兆通り あるみたいです。

ご質問等ございましたら、hamano@fusic.co.jpまで、ご連絡ください。

hamano

hamano

I'm a software engineer in Fukuoka, Japan. Recently, I am developing machine learning models using TensorFlow, and also developing Web services by using PHP.