hamano 準備
準備として、各種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
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.