시뮬레이티드 어닐링 (Simulated Annealing)

시뮬레이티드 어닐링 (SA, simulated annealing)은 1983년에 커크패트릭(Kirkpatrick)이 개발한 최적화 알고리즘으로서 물리적 어닐링 과정을 최적화의 관점에서 모방하여 알고리즘을 개발했다고 한다. 시뮬레이티드 어닐링(SA)은 처음에 순회외판원문제(TSP, traveling salesman problem)와 같은 조합 최적화 문제를 해결하기 위해 도입되었지만 최근에는 연속형 최적화 문제 해결에도 적용되고 있다.

물리적 어닐링은 금속을 가열하여 액체 상태로 만든 다음 규칙적인 결정성 구조를 형성할 때까지 천천히 냉각시키는 방법이다. 냉각 과정에는 냉각이 시작되는 초기 온도와 냉각 속도 등 두가지 주요 파라미터가 있는데, 이 값에 따라 고체 상태의 금속의 강성과 원자 구조의 규칙성이 결정된다. 초기 온도가 충분히 높지 않으면 액체 상태에서 금속의 원자가 최소 에너지 구조로 위치를 재배열할 수 있는 자유도가 충분하지 않고, 냉각 시간이 충분하게 주어지지 않으면 금속의 원자가 최소의 내부 에너지로 규칙적으로 정렬하여 결정 격자를 형성 할 수 있는 시간이 부족해진다.

 

 

'시스템'이 특정 상태에 있을 확률을 해당 상태의 에너지와 시스템의 온도의 함수로 나타내는 확률 분포를 볼츠만(Boltzman) 분포라고 하며 다음과 같이 주어진다.

 

$$\begin{array}{}\text{(1)}& p_i\propto \exp (-\frac{E_i}{k_{B}T})\end{array}$$

 

여기서 시스템이란 충분한 수의 원자의 집합을 의미하며, $p_i$ 는 시스템이 상태 $i$ 에 있을 확률, $E_i$ 는 해당 상태의 에너지, $k_B$ 는 볼츠만 상수, $T$ 는 온도다. 이 분포에 의하면 에너지가 낮은 상태의 확률이 높고 $T$ 가 낮을 수록 낮은 에너지 상태에 확률이 집중되어 있다. 상태의 에너지 $E_i$ 를 양자화된 값으로 보지않고 연속적인 값으로 가정하면 식 (1)은 확률밀도함수(pdf)로 해석해도 무방하다.

 

 

SA 알고리즘은 물리적 어닐링을 시뮬레이션하는 알고리즘으로서 어닐링이 진행되는 금속의 물리적 구조의 변화를 재현함으로서 최소의 에너지를 갖는 해(solution)를 구하는 것을 목표로 한다.

다음과 같이 제약조건이 없는 일반적인 최적화 문제를 보자.

 

$$\begin{array}{}\text{(2)}& {\mathbf{x}}^{\star }=\arg \underset{\mathbf{x}}{min}f(\mathbf{x})\end{array}$$

 

여기서 $\mathbf{x}\in {\mathbb{R}}^n$ 은 최적화 변수이고, $f(\mathbf{x})$ 는 목적함수(objective function)이다. 이와 같은 최적화 문제에 대한 SA 알고리즘은 다음과 같다.

0. 온도(temperature)의 초기값 $T_1$ 을 정한다. 그리고 최적화 문제의 첫번째 해(solution) ${\mathbf{x}}^{(1)}$ 을 임의로 선정한다. $E_1=f({\mathbf{x}}^{(1)})$

다음을 $T_i=T_{min}$ 이 될 때까지 반복한다.

1. 후보(candidate) 해를 추출한다.
이에 대해서는 다양한 방법이 제안되었는데 대개 현재 해 ${\mathbf{x}}^{(i)}$ 를 조건으로 하여 일정 범위 내에서 무작위로 선정하든가 아니면 대칭인 제안(proposal) 확률 분포에서 추출한다.

              $\mathbf{x}\sim \text{generate }({\mathbf{x}}^{(i)})$
        $\text{or }\mathbf{x}\sim q(\mathbf{x}|{\mathbf{x}}^{(i)})$

2. 만약 $\mathrm{\Delta }E=E-E_i=f(\mathbf{x})-f({\mathbf{x}}^{(i)})<0$ 이면 후보 해를 새로운 해로 수락한다 (${\mathbf{x}}^{(i+1)}=\mathbf{x}$ ).

만약 $\mathrm{\Delta }E\ge 0$ 이면,

     2-1. 후보 해의 수락 확률(acceptance probability)을 계산한다.

 

$$\begin{array}{rl}\alpha & =min(1,\frac{\exp (-\frac{E}{T_i})}{\exp (-\frac{E_i}{T_i})})\\ \\ & =min(1,\exp (-\frac{E-E_i}{T_i}))\\ \\ & =min(1,\exp (-\frac{\mathrm{\Delta }E}{T_i}))\end{array}$$

 

     2-2. 확률 $\alpha$ 로 후보 해 $\mathbf{x}$ 를 새로운 해로 수락하거나 거절한다.
즉 균등분포 $\mathcal{U}[0,1]$ 에서 샘플 $u$ 를 추출한 후, $u\sim \mathcal{U}[0,1]$, 만약 $u<\alpha$ 이면 후보 해를 새로운 해로 수락하고 아니면 거절한다.

3. 선택한 냉각(cooling) 스케줄에 따라 다음 온도 $T_{i+1}$ 를 낮춘다.

위 알고리즘에서 보듯이 SA는 해를 수정하는 과정과 해를 선택적으로 수용하는 수락 확률로 구성되는 메트로폴리스(Metropolis) 기법을 기반으로 한다. 메트로폴리스 알고리즘에서는 제안 확률밀도함수가 대칭이라고 가정하기 때문에 수락 확률이 목표 확률밀도함수만의 함수가 된다. 여기서 목표 확률밀도함수 또는 확률 분포에 해당하는 것은 식 (1)의 볼츠만 분포다. 또한 에너지는 $E=f(\mathbf{x})$ 로서 목적함수의 값에 해당한다.

SA는 물리적 어닐링 과정을 모방하는 데서 유래하여 볼츠만 분포를 목표 확률 분포로 간주하지만 일반적인 확률 분포 $p(\mathbf{x})=\frac{\tilde{p}(\mathbf{x}))}{Z}$ 를 목표 확률 분포로 정해도 된다. 이 경우 수락 확률을 다음과 같이 계산하면 된다.

 

$$\begin{array}{r}\alpha =min(1,\frac{{\tilde{p}}^{\frac{1}{T_i}}(\mathbf{x})}{{\tilde{p}}^{\frac{1}{T_i}}({\mathbf{x}}^{(i)})})\end{array}$$

 

여기서 ${\tilde{p}}^{\frac{1}{T_i}}(\mathbf{x})$ 는 $T_i\to 0$ 일 수록 $\tilde{p}(\mathbf{x})$ 의 최대값을 더욱 과장시키는 확률밀도함수다.

 

 

SA 알고리즘과 메트로폴리스 알고리즘(https://pasus.tistory.com/367)을 비교해 보면, SA 알고리즘이 목표 확률밀도함수 $\tilde{p}(\mathbf{x})$ 대신에 ${\tilde{p}}^{\frac{1}{T_i}}(\mathbf{x})$ 을 사용한 것을 제외하고는 두 알고리즘이 거의 같다는 것을 알 수 있다.

메트로폴리스 알고리즘에 의하면 이론적으로는 특정 온도 $T_i$ 에서 마르코프 체인이 정상(stationary) 상태에 도달할 때까지 기다려야 하나 SA 알고리즘에서는 온도 $T_i$ 를 서서히 줄임으로써 준평형 상태에 도달할 수 있도록 한다.

SA 알고리즘에서 온도 $T_i$ 를 처음부터 매우 작게 설정하면 되지 않을까 하는 의문이 들 수도 있는 데, 이 경우에는 함수의 최대값 또는 최소값 근처에서 샘플링될 확률이 극도가 작아지기 때문에 효율이 매우 떨어진다. 따라서 서서히 온도를 감소시키면서 해당 지역에서 샘플링될 확률을 점차로 높이는 것이다. 이는 물리적인 어닐링 과정과 유사하다.

 

 

냉각 또는 어닐링 스케줄, 즉 온도를 낮추는 방법은 기존 온도에서 일정한 온도를 빼주는 방법 ($T_{i+1}=T_i-T_d$)과 기존 온도에 스케일 팩터를 곱하는 방법 ($T_{i+1}=R_d{T}_i$)이 있다 보통 스케일 팩터를 사용하는 방법이 성능이 좋다고 알려져 있다.

효율적인 어닐링 알고리즘을 얻으려면 적절한 제안 분포와 적절한 냉각 스케줄을 선택하는 것이 매우 중요하다. 문헌에 보고된 부정적인 SA 결과의 대부분은 잘못된 제안 분포 설계에서 비롯된 경우가 많다고 한다.

 

 

다음은 PSO에 관한 게시글(https://pasus.tistory.com/357)에서 예제로 풀었던 최적화 문제다. 이번에는 SA로 풀어보자.

 

$$\begin{array}{rl}& {\mathbf{x}}^{\star }=\arg \underset{\mathbf{x}}{min}f(\mathbf{x})\\ & f(\mathbf{x})=(x_1-1{)}^2+(x_2-2{)}^2+(x_3+3{)}^2,\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]\end{array}$$

 

답이 ${\mathbf{x}}^{\star }=[1,2,-3{]}^T$ 로 정해져 있으니 SA 결과와 비교하면 좋을 것이다. 초기 온도를 $T_1=1$, 냉각 스케줄을 $T_{i+1}=0.995T_i$, 제안 확률밀도함수를 표준편차가 $1$ 인 가우시안으로 했을 때 실행 결과는 다음과 같다.

 

iteration 1 : cost = 33.52190 
iteration 2 : cost = 33.52190 
iteration 3 : cost = 33.52190 
iteration 4 : cost = 33.52190 
iteration 5 : cost = 21.72813 
iteration 6 : cost = 21.72813 
iteration 7 : cost = 21.72813
...
iteration 998 : cost = 0.00782 
iteration 999 : cost = 0.00782 
iteration 1000 : cost = 0.00782

solution is 
  0.97813
  2.08568
 -2.99964

 

 

매트랩 코드는 다음과 같다.

 

%
% Simulated Annealing for argmin f(x)
%
% (c) st.watermelon

clear;

% Gaussian proposal pdf with mu-mean and 1-std
q_pdf = @(x, mu) 1/sqrt(2*pi)*exp(-(x-mu)^2/2);

% parameters
T0 = 1;       % initial temperature
T = T0;

cooling_rate = 0.995;
max_iter = 1000;

% 0. initialize
x = randn(3,1);
fn = objective(x);

for t = 1 : max_iter

    fn_hist(t) = fn';

    % 1. sampling candidate from proposal q(x_new|x_t)
    x_new = randn(3,1) + x;
    
    % 2. accept or reject
    fn_new = objective(x_new);
    dE = fn_new - fn; % to downhill 
    if dE < 0
        x = x_new;
        fn = fn_new;
    else
        % 2-1. acceptance probability
        alpha = min(1, exp(-dE/T));
    
        % 2-2. accept or reject
        u = rand(1);
        if u < alpha
            x = x_new;
            fn = fn_new;
        end
    end

    % 3. update temperature
    T = T * cooling_rate;

    if T < 1e-6
        break;
    end

    fprintf('iteration %d : cost = %.5f \n', t, fn);

end


% result
semilogy(fn_hist);
fprintf('\n solution is \n %8.5f\n %8.5f\n %8.5f\n\n', x);

%%
function out = objective(x)

out = (x(1)-1).^2 + (x(2)-2).^2 + (x(3)+3).^2;

end