[수학적 최적화] 시뮬레이션 담금질 Simulated Annealing을 활용한 외판원 문제 Travelling salesman problem 문제 풀기

지난 포스팅에서 생성한 시뮬레이션 담금질 프로그램을 활용하여, 대표적인 조합 최적화combinatorial optimization 문제인 외판원 문제travelling salesman problem을 풀어 보도록 하겠습니다. 시뮬레이션 담금질에 대해서 설명한 지난 포스팅을 아직 읽지 않으신 분이라면, 이전 포스팅을 먼저 읽고 이 포스팅을 읽으시길 바랍니다. 

https://studyingrabbit.tistory.com/100

[수학적 최적화] 시뮬레이션 담금질 Simulated Annealing의 이해와 파이썬 코드

외판원 문제

 

외판원 문제에 대해서는 한 번 쯤은 들어았을 것이라 생각합니다. 외판원 문제는

지도 상에 여러 도시들이 있고, 한 도시에서 다른 도시 사이의 거리가 주어져 있을 때, 모든 도시들을 단 한 번만 방문하고 여행을 시작 했던 도시로 돌아올 때, 이동 거리를 최소로 하는 도시 방문 순서를 찾는 문제

입니다. 물론 문제의 성격에 따라서 약간의 차이가 있을 수 있습니다. 

 

외판원 문제는 일상생활에서도 겪는 문제입니다. 특히 택배 기사분들은 매일 겪는 문제인데요, 배달의 목적지가 여러곳인 경우, 어떤 순서로 배송을 해야 가장 빠르게 물건을 배송할 수 있을지 고민하실 것 입니다. (약간은 다른 문제이긴 하지만) 여러 마을을 지나는 버스나 셔틀의 경로를 결정할 때, 외판원 문제를 풀게 됩니다. 

 

문제의 정확한 정의

 

지도상에 $N$개의 도시의 위치와 임의의 두 도시간 거리가 주어졌다고 가정하겠습니다. 나중에 수학적 최적화 문제를 풀기 위해서 명확하게 정의를 하면, $1$번 부터 $N$번 까지의 도시를 $x_i, i=1, 2, 3, ...N$ 이라고 표현하고, $i$ 번째 도시와 $j$번째 도시 사이의 거리를 $d_{ij}$ 라고 하겠습니다. 또한, $N$개의 도시를 방문 하는 경로를 수열 $a_n$으로 쓰도록 하겠습니다. $a_n = (a_1, a_2, ..., a_N)$ 이며, $a_i \in \{ 1, 2, 3, .., N\}$ 는 $i$번째 방문하는 도시 입니다. 또한 편의상, $a_{N+1} = a_1$로 정의하겠습니다. 왜냐하면, 모든 도시를 방문한 이후에, 처음에 출발했던 도시로 돌아가야 하기 때문입니다. 예를 들어 $N=4$일 때, 경로 $a_n = (2, 1, 3, 4)$는 $x_2 \rightarrow x_1 \rightarrow x_3 \rightarrow x_4$ 순서로 도시를 방문하여 최종적으로 $x_2$로 돌아오는 경로를 의미합니다. 임의의 경로 $a_n$가 주어졌을 때, 이 경로에 대응하는 거리는 
$$D = D(a_n) = \sum_{i = 1}^{N} d_{a_i, a_{i+1}}$$

로 주어집니다. 위에서 사용한 예시를 사용하면,

$$D = D(2, 1, 3, 4) = d_{2,1} + d_{1,3} + d_{3,4} + d_{4,2}$$

이 됩니다. 

 

외판원 문제는 매우 풀기 어려운 문제다

 

$N$개의 도시가 주어졌을 때, 이 도시를 방문하는 방법은 $N!$ 개가 됩니다. 여기서 $!$는 팩토리얼을 의미 합니다. 단순하게 생각했을 때, 시작 하는 도시를 선택했다면, 그 다음에 방문하는 도시는 시작 하는 도시를 제외한 $N-1$개의 도시가 됩니다. 그 다음에 방문할 수 있는 도시는 $N-2$개가 되고.... 이렇게 모든 가능한 조합을 생각하면 $N!$ 개가 됩니다. 물론 시작하는 위치로 다시 돌아 오는 경우를 생각하기 때문에, 경로를 의미하는 수열 $a_n$이 서로 다르다고 하더라도, 실제로 동일한 경로를 의미하는 경우도 있습니다. 예를들어, $a = (1, 2, 3, 4)$ 와 $b = (2, 3, 4, 1)$는 동일한 경로를 의미합니다. 따라서 실제 가능한 경로의 갯수는 $N!$ 보다 작습니다. 

 

$N$이 커짐에 따라서 팩토리얼 $N!$는 매우 크게 증가하는 함수 입니다. 만일 10개의 도시가 주어졌다면, 최적의 경로를 찾는 문제는 총 $10! = 3,528,800$ 개 입니다. 도시가 10개 밖에 되지 않는데도 최적의 경로의 후보가 되는 경로는 총 300만개가 넘습니다! 300만개의 경로에 대해서 경로에 대응하는 거리를 계산하는 것은 무척이나 어렵습니다. "300만개 정도는 계산해 볼 수 있겠다" 라는 생각이 들 수 도 있습니다. 그렇다면 10개의 도시가 더 추가되어, 총 20개의 도시가 있다면 어떨까요? 이 때 모든 경로의 갯수는$20! = 2.4 \times 10^{18}$ 개가 됩니다! 이건 계산할 수 있는 숫자가 아닙니다! 따라서 단순히 "모든 경로에 대한 거리를 계산 한 후, 이 중에서 최소값을 선택하는 방법"은 통하지 않습니다.

 

외판원 문제는 계산 복잡성 이론에 따르면 NP-hard에 속하는 문제 입니다. 이 문제를 풀기 위해 다양한 방법이 시도되고 있습니다. 우리는 이번 포스팅에서 이 문제를 이전 포스팅에서 공부한 바 있는 시뮬레이션 담금질을 통해서 풀어 보고자 합니다. 시뮬레이션 담금질을 통해서 얻어진 해가 실제 최적의 해라는 것을 보장할 수는 없지만 현실적으로 볼 때 "꽤 괜찮은 해" 정도는 됩니다. 

 

시뮬레이션 담금질 방법을 통한 외판원 문제 풀기

 

1. 목적 함수의 정의역

 

지난 포스팅에서 개발한 시뮬레이션 담금질 코드를 약간만 수정하면, 외판원 문제에 적용할 수 있습니다. 지난 포스팅에서 다룬 문제와 이번 포스팅에서 다루는 외판원 문제의 차이는 목적 함수의 정의역이 다르다는 것 입니다. 그 밖에 차이는 없습니다. 물론 주어진 목적 함수의 정의역의 한 점 (상태) 에 대응하는 목적 함수를 계산하는 방식이 조금 다르긴 하지만, 이는 큰 문제는 되지 않습니다. 그냥 목적 함수를 계산하는 부분의 코딩을 각각에 맞게 하면 됩니다. 

 

목적 함수의 정의역은 2차원 유클리드 공간의 부분집합이었습니다. 즉 $(x, y)$로 기술이 됩니다. 이번 문제의 목적 함수의 정의역은 경로입니다. $N$개의 도시가 주어진 경우, 정의역의 한 상태는 앞에서 기술한 대로 $a_n = (a_1, a_2, ...,. a_N)$으로 기술 됩니다. 여기서 $a_i \in \{ 1, 2, 3, ..., N\}$ 입니다. 목적 함수의 정의역이 다르니, 이 정의역에 속한 한 상태에서 다른 상태로의 "이동"을 정의하는 방법 역시 바뀌어야 합니다. 이전 포스팅의 문제에서는

$$(x_{old}, y_{old}) \rightarrow (x_{new}, y_{new}) = (x_{old} + 0.1 \cdot r_1, x_{old} + 0.1 \cdot r_2), r_1, r_2 \in [-1, 1]$$

와 같이 한 상태 $(x_{old}, y_{old})$ 에서 새로운 상태 $(x_{new}, y_{new})$를 정의하였습니다. 이는 현재의 상태에서 조금 변화를 주어 새로운 상태를 얻은 것 입니다. 비슷한 개념을 외판원 문제에도 적용할 수 있습니다. 현재의 경로 $a_n$이 주어졌을 때, 이 경로를 약간 변화 시키는 새로운 경로 $b_n$을 만들어 내면 되는데, 생각해 볼 수 있는 간단한 방법은, "임의의 두 도시(수열상의 숫자)를 선택하고, 이 두 도시의 순서를 서로 바꾼는 것" 입니다. 예를들어 $a_n = (a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8)$가 주어져 있을 때, 1부터 8사이 숫자 중 두 개를 고릅니다. 예를들어 $(4, 7)$이 선택 되었다고 하겠습니다. 그러면 $a_4, a_7$의 순서를 바꾸는 것 입니다. 즉 $b_n = (a_1, a_2, a_3, a_7, a_5, a_6, a_4, a_8)$이 됩니다.

 

조금 더 명확하게 수학적으로 표시하자면, $1$ 부터 $N$ 까지의 숫자 중, 서로 다른 두 숫자를 임의의로 선택하여 $p, q$라고 했을 때,$b_n = (b_1, b_2, .., b_N) $에 대해서,

$$\begin{equation}
      b_i=\begin{cases}
        a_q, \text{ if } i=p\\
        a_p, \text{ if } i=q \\
        a_i, \text{ otherwise }
      \end{cases}
    \end{equation}$$

가 됩니다. 

 

2. 목적 함수

 

외판원 문제의 목적 함수는 주어진 경로에 대응하는 거리를 구하는 것 입니다. 임이의 두 도시 쌍 $(i,j)$에 대해서 이 두 도시간의 거리 $d_{i,j}$가 주어진다면, 위에서 설명한 식에 따라서 거리를 구할 수 있습니다. 우리가 풀어 볼 문제에서는 각 도시가 2차원 공간상에 주어졌을 때, 유클리드 거리를 계산하여 $d_{i,j}$를 정의할 것 입니다. 즉, $i, j$ 번째 도시의 위치를 $(x_i, y_i), (x_j, y_j)$이라고 한다면, $d_{i,j} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}$ 입니다. 

 

3. 문제 세팅

 

이번 포스팅에서 실제로 풀어볼 문제에 대해 설명을 하도록 하겠습니다. 

 

$30$개의 "도시"가 위 그림의 녹색 점과 같이 원주상의 한 점에 놓여 있다고 하겠습니다. 맨 오른쪽에 위치한 점을 $x_1$ 이라고 하고 반 시계 방향대로 $x_2, x_3, ..., x_30$ 이라고 하도록 하겠습니다. "지나치게 이상한 배치" 라고 하실 수 있지만, 사실 도시의 위치 $(x, y)$가 중요한 것이 아니라, 두 도시간의 거리 $d_{i,j}$가 중요한 요소 입니다. 위와 같이 도시의 위치를 배치한 것은 위와 같이 도시를 배치하면 최적의 경로(거리가 가장 짧은 경로)를 쉽게 알 수 있기 때문입니다. 위와 같이 도시가 배치되어 있다면, 최적의 경로는 누가 보더라도 $a = (1, 2, 3, 4, ..., 27, 28, 29, 30)$ 입니다. 만일 위와 같이 배치 하지 않고 임의의 30개의 점을 배치한다면, 어떤 경로가 최적의 경로인지를 알 수 없고, 따라서 우리가 시뮬레이션 담금질로 부터 구한 해가 실제로 최적의 경로인지를 확인할 방법이 없습니다. 

 

$N = 30$의 외판원 문제이니, 가능한 총 경로의 갯수는 $30! = 2.6 \times 10^{32}$ 개나 됩니다! 시뮬레이션 담금질이 이 많은 경로 가운데서 최적의 경로를 찾아 낼 수 있을지 궁금합니다. 

 

4. 풀어 보기

 

위 문제에 대한 파이썬 코드 입니다. 

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation
import random

plt.style.use('dark_background')

def objective_tsp(x, P):
    dist = 0.0
    for index in range(len(x)):
        dx = P[x[index]] - P[x[(index+1) % len(P)]]
        d = np.linalg.norm(dx)
        dist += d
    return dist


def transition_probability(e, e_new, T):
    de = e_new - e
    p = np.exp(-de/T)
    return p

def transition_or_not(p):
    r = np.random.random()
    if p > r:
        return True
    else:
        return False

def update_tsp(num, list_line, list_E, list_T, line, time_text):
    line.set_data(list_line[num][:,0], list_line[num][:,1])
    time_text.set_text("Iteration = {}\nObjective = {:.4f}\nTemperature = {:.4f}".format(num, list_E[num], list_T[num]))
def tsp():

    Pt = 30
    list_theta = np.linspace(0, 2*np.pi, Pt)[:Pt-1]
    P = [[0.5+0.4*np.cos(theta), 0.5+0.4*np.sin(theta)] for theta in list_theta]
    P = np.array(P)
    E_opt = objective_tsp(range(Pt-1), P)

    P = P[:]

    T = 10
    Niteration = 10000
    T_reduce = 10 ** (-3 / Niteration)
    N_skip = Niteration//400

    X = list(range(len(P)))
    np.random.shuffle(X)

    list_X, list_E, list_T = [], [], []

    for iteration in range(Niteration):
        X_new = list(X)

        a, b = random.sample(range(len(X)), 2)
        X_new[a], X_new[b] = X_new[b], X_new[a]

        e = objective_tsp(X, P)
        e_new = objective_tsp(X_new, P)
        p = transition_probability(e, e_new, T)
        bool_transition = transition_or_not(p)

        if bool_transition:
            X = X_new
        T *= T_reduce

        list_X.append(X)
        list_E.append(e)
        list_T.append(T)

        print(iteration, e, T, E_opt)

    list_X = np.array(list_X)

    list_line = []
    for X in list_X:
        line = []
        for index in range(len(X)):
            p = [P[X[index]][0], P[X[index]][1]]
            line.append(p)
        line.append([P[X[0]][0], P[X[0]][1]])
        list_line.append(line)

    list_line = np.array(list_line)

    fig = plt.figure(figsize=(8, 8))
    ax = plt.axes()

    point = ax.scatter(P[:,0], P[:,1], s=50, color='lime', zorder=1)
    line, = ax.plot([], [], color='lime', lw=2)
    time_text = ax.text(0.65, 0.05, '', transform=ax.transAxes, fontsize='large')
    plt.xlabel('X')
    plt.ylabel('Y')

    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)

    list_line = list_line[::N_skip]
    list_E = list_E[::N_skip]
    list_T = list_T[::N_skip]

    ani = animation.FuncAnimation(fig, update_tsp, len(list_E), fargs=(list_line, list_E, list_T, line, time_text),
                                  interval=10, blit=False)
    plt.show()
    plt.close('all')

    plt.plot(range(len(list_E)), list_E, color = 'yellow', label = 'SA')
    plt.axhline(E_opt, color = 'lime', lw = 1, ls = '--', label = 'exact')
    plt.xlabel("Iteration")
    plt.ylabel("Objective")
    plt.legend()
    plt.grid(lw = 1, color = 'gray')
    plt.show()

if __name__ == '__main__':
    tsp()

코드에 대한 설명과, 지난 포스팅 대비 이번 포스팅에서 달라진 부분에 대해서는 앞에서 설명을 했기 때문에 코드 부분 부분에 대한 설명은 생략 하도록 하겠습니다. 

 

위 코드를 실행하면, 아래와 같은 결과를 얻게 됩니다. 물론 프로그램에는 임의성(randomness)를 포함하고 있기 때문에, 정확히 똑같은 결과가 얻어지는 것은 아닙니다. 

https://youtu.be/4gLzaDq07yg

반복의 초기에는 녹색 실선들이 어지럽게 펼쳐져 있다가 끝 부분이 되면 가장 가까이에 있는 점을 잇는 선분 형태로 최적의 경로를 찾아 갑니다. 영상의 마지막에는 반복 시행에 따른 목적 함수 (해당 경로의 거리)가 나옵니다. 이 시뮬레이션은 총 10,000번의 반복을 시행한 것이며, 동영상에서는 시각화의 편의를 위해서 총 400개의 스냅샷만 표현 됩니다. 즉, 25 번의 반복 중에서 단 한 번의 이미지만 보여주는 것 입니다. 

 

반복의 시행에 따라 찾아진 경로가 어떻게 최적의 경로로 옮아 가는지 확인해 보도록 하겠습니다. 

 

시뮬레이션의 초반 입니다. 시작 경로는 $30$개의 도시를 완전 무작위로 순서를 매긴 것이기 때문에, 위와 같이 경로를 잇는 선분이 중구난방으로 표현 됩니다. 당연히 경로의 거리는 깁니다. 

 

반복 시행을 하면 할 수록 최적의 경로를 찾아 갑니다. 위 스냅샷에서는 대략 한 절반 조금 못 되는 정도로 최적의 경로(바로 옆에 있는 도시와 연결 됨)를 찾을 것 같고, 절반 정도는 최적의 경로를 찾지 못 했습니다. 길게 표시 되는 선분이 5개 정도 있는 것을 봐서는, 아직까지는 정답에 비슷한 경로를 찾지 못하는 경우도 좀 있습니다. 

 

반복 시행을 조금씩 더 해가면 최적의 경로와 매우 유사한 경로를 찾게 됩니다. 위 경우에는 아직 5곳에서는 최적의 경로를 찾지 못 했습니다. 목적 함수의 변화량이 대략 1 정도의 수준이고, 온도 값이 대략 0.01 오더 이니, 이 정도 시점부터는 국소 최소점을 벗어나려는 시도 보다는 주어진 최소점 근처에서 보다 작은 목적 함수 값을 찾으려는 시도를 합니다. 

 

동영상 상에서 약 300번째 반복, 전체 반복 횟수에서는 약 7,500번째 반복에서 최적 해를 찾게 됩니다. 이후에도 새로운 경로를 찾으려는 시도를 당연히 하지만, 목적 함수의 값이 증가하는 시도이기에 새로운 경로로 이동하지는 않습니다. 

 

반복 시행에 따른 목적 함수의 변화를 보여주는 그래프 입니다. 반복 시행 초기에는 온도가 높기 때문에 이런 저런 상태를 막 오고 갑니다. 따라서 목적 함수의 값 역시 중구난방으로 튀게 됩니다. 그러다가 반복 시행이 늘어날 수록 목적 함수의 값은 극적으로 낮아집니다. 이 과정에서도 목적 함수가 증가하는 시행이 있긴 하지만, 그 횟수는 반복의 초기 보다는 작습니다. 

 

제가 시행했을 때는 위와 같이 성공적으로 최적의 경로를 찾았습니다. 그러나 앞서 설명한대로 코드에는 무작위성이 있기 때문에 항상 최적 경로를 찾을 수 있는 것은 아닙니다. 이 경우에는 반복 시행의 횟수를 늘리거나, 온도를 좀 더 키우는 방식으로 최적 경로를 찾을 수 있게 할 수 있습니다. 

 

시뮬레이션 담금질이 외판원 문제를 제대로 풀었다!

 

30개의 도시에 대한 외판원 문제의 경로는 총 대략 $10^{32}$ 개나 됩니다. 이렇게 많은 최적 경로의 후보 중에서, 우리의 시뮬레이션 담금질은 최적의 경로를 찾아 냈습니다. 그러나 어찌 생각하면 $N = 30$은 꾀나 쉬운 문제라고 생각할 수 있습니다. 도시의 갯수를 키웠을 때도, 우리의 코드가 제대로된 값을 얻어 낼 수 있는지는 각자 테스트 해 보면 재미있을 것 같습니다.