巡回セールスマン問題

総移動距離の計算

次に、総移動距離の計算をする関数を用意します。これが訪問順の評価指標となります。距離を何度も計算するのはもったいないので、先に都市間の距離を計算しておきます。

x = city_xy[:, 0]
y = city_xy[:, 1]
distance_matrix = np.sqrt((x[:, np.newaxis] - x[np.newaxis, :]) ** 2 +
                          (y[:, np.newaxis] - y[np.newaxis, :]) ** 2)

# 見やすくするために10都市分だけ出力。
# distance_matrix[i, j]が都市iと都市jの距離。
print(np.round(distance_matrix[:10, :10]))

[[   0.   74.   34.   94.   61.   94.   91.   59.   28.   71.]
 [  74.    0.   54.   44.   81.   21.   67.   22.   47.   29.]
 [  34.   54.    0.   61.   36.   75.   57.   50.   23.   64.]
 [  94.   44.   61.    0.   67.   52.   32.   63.   70.   72.]
 [  61.   81.   36.   67.    0.  101.   45.   83.   59.   98.]
 [  94.   21.   75.   52.  101.    0.   81.   36.   66.   33.]
 [  91.   67.   57.   32.   45.   81.    0.   81.   75.   94.]
 [  59.   22.   50.   63.   83.   36.   81.    0.   33.   15.]
 [  28.   47.   23.   70.   59.   66.   75.   33.    0.   47.]
 [  71.   29.   64.   72.   98.   33.   94.   15.   47.    0.]]

def calculate_total_distance(order, distance_matrix):
    """Calculate total distance traveled for given visit order"""
    idx_from = np.array(order)
    idx_to = np.array(order[1:] + [order[0]])
    distance_arr = distance_matrix[idx_from, idx_to]

    return np.sum(distance_arr)

# 試しに距離を計算してみる
test_order = list(np.random.permutation(N))
print('訪問順序 = {}'.format(test_order))

total = calculate_total_distance(test_order, distance_matrix)
print('総移動距離 = {}'.format(total))

訪問順序 = [9, 14, 0, 12, 15, 6, 1, 11, 16, 2, 3, 7, 4, 18, 19, 5, 17, 10, 13, 8]
総移動距離 = 1079.4019034407843

訪問経路の可視化

経路を可視化するコードも先に準備しておきます。こういう計算をする時に、数字だけ見て作業をするよりも色々気付きがあるので、可視化大事です。あと、「やってる感」が出てやる気が出ます。（笑）

def visualize_visit_order(order, city_xy):
    """Visualize traveling path for given visit order"""
    route = np.array(order + [order[0]])  # add point of departure
    x_arr = city_xy[:, 0][route]
    y_arr = city_xy[:, 1][route]

    plt.figure(figsize=(4, 4))
    plt.plot(x_arr, y_arr, 'o-')
    plt.show()

# 可視化を試してみる
visualize_visit_order(test_order, city_xy)

2-opt近傍

この「現在の解に少し変化を加える」方法は問題毎にさまざまな方法が提案されています。ここでは、循環セールスマン問題でよく使われる近傍の例として、2-opt近傍を使います。2-opt近傍は、経路中の2本の移動を交換するものです。素直に実装すると、交換前後でそれぞれ総移動距離を計算して比較すれば良いのですが、結局交換した部分のみ考えても同じことなので、交換による総移動距離の差分を計算する関数を準備します。

def calculate_2opt_exchange_cost(visit_order, i, j, distance_matrix):
    """Calculate the difference of cost by applying given 2-opt exchange"""
    n_cities = len(visit_order)
    a, b = visit_order[i], visit_order[(i + 1) % n_cities]
    c, d = visit_order[j], visit_order[(j + 1) % n_cities]

    cost_before = distance_matrix[a, b] + distance_matrix[c, d]
    cost_after = distance_matrix[a, c] + distance_matrix[b, d]
    return cost_after - cost_before

以下が、２つのパスの交換後の訪問順序を計算するコードになります。

def apply_2opt_exchange(visit_order, i, j):
    """Apply 2-opt exhanging on visit order"""

    tmp = visit_order[i + 1: j + 1]
    tmp.reverse()
    visit_order[i + 1: j + 1] = tmp

    return visit_order

近傍探索の実装

現状の訪問経路の、各２パスを入れ替える操作を全組み合わせで実施します。全通りを計算しておき、もっとも総移動距離を減らせる交換を実際に適用します。これ以上改善できなければNoneを返すことにします。

def improve_with_2opt(visit_order, distance_matrix):
    """Check all 2-opt neighbors and improve the visit order"""
    n_cities = len(visit_order)
    cost_diff_best = 0.0
    i_best, j_best = None, None

    for i in range(0, n_cities - 2):
        for j in range(i + 2, n_cities):
            if i == 0 and j == n_cities - 1:
                continue

            cost_diff = calculate_2opt_exchange_cost(
                visit_order, i, j, distance_matrix)

            if cost_diff < cost_diff_best:
                cost_diff_best = cost_diff
                i_best, j_best = i, j

    if cost_diff_best < 0.0:
        visit_order_new = apply_2opt_exchange(visit_order, i_best, j_best)
        return visit_order_new
    else:
        return None

上で実装した2-optによる近傍探索を、改善が見られる限り繰り返します。

def local_search(visit_order, distance_matrix, improve_func):
    """Main procedure of local search"""
    cost_total = calculate_total_distance(visit_order, distance_matrix)

    while True:
        improved = improve_func(visit_order, distance_matrix)
        if not improved:
            break

        visit_order = improved

    return visit_order

これで近傍探索のコードは完成です。早速、近傍探索を動かしてみましょう。

# 適当に初期解を生成
test_order = list(np.random.permutation(N))
visualize_visit_order(test_order, city_xy)
total_distance = calculate_total_distance(test_order, distance_matrix)
print('初期解の総移動距離 = {}'.format(total_distance))

# 近傍を計算
improved = local_search(test_order, distance_matrix, improve_with_2opt)
visualize_visit_order(improved, city_xy)
total_distance = calculate_total_distance(improved, distance_matrix)
print('近傍探索適用後の総移動距離 = {}'.format(total_distance))

初期解の総移動距離 = 1089.0058941597254

近傍探索適用後の総移動距離 = 420.08788455293126

ランダムの初期解に従って各都市を回ってこい！と言われたら心が折れそうですよね。それに対して近傍探索で改善を行った経路の方は、近い都市を順番にまわるような経路に改善されているのが分かります。