模拟退火算法原理

• 模拟退火算法是一种通用概率算法，用来在一个大的搜寻空间内寻找问题最优解。
• 其思想源于固体的退火过程：将固体加热至足够高的温度，缓慢冷却，其内能就由很大缓慢趋于内能最小。
• 足够高的温度对应随机的解，缓慢冷却对应对解进行一次随机扰动，内能对应目标函数。即将随机解不断进行扰动，根据目标函数变化以一定概率接受解，不断重复，解就会趋近于最优解。

• 接受新解概率——Metropolis准则：

  • 设从当前状态 $i$ 生成新状态 $j$ ，若 $E_j<E_i$，则新状态 $j$ 作为当前状态，否则，以概率 $exp[\frac{-(E_j-E_i)}{k×t}]$ 接受当前状态。 $k$ 为常数，通常取1； $t$ 是此时温度。
  • 内能对应由于目标函数（求最小值，如果求最大值可加负号求最小），如果新解优于旧解，则接受，否则以一定概率接受。
  • 可以看出，概率和新旧解的差值有关，相差越小，越可能被接受，即有概率跳出此时所在的局部，以便在全局寻找。

模拟退火算法步骤

符号说明

算法步骤

• 令T = T₀，随机生成初始解 $x_0$，计算目标函数 $E(x_0)$
• 根据实际情况对当前解进行扰动，产生新解，再次计算目标函数。
• 比较目标函数，根据Metropolis准则接受新解。
• 重复扰动 L_k 次，达到次数后再按系数α降低温度。
• 判断T是否达到T_f ，达到则终止算法，否则继续扰动。

注： 一般会把退火过程中碰到的最优解也保存下来，与最后一个解一同输出，因为有可能在以一定概率接受较差解的时候接受了较差解。

扰动方法

• 对应实际问题，都有不同的扰动方法以及不同的目标函数设计方法。
• 解x一般是能够表示一种实际方法的量，比如一个数组，根据x可以得出目标函数E的值。
• 根据实际情况如何变动，设计对解x进行改变，也就实现了扰动。
• 产生的随机扰动要保证可能取遍所有可能情况。

实例——旅行商问题

问题描述

一名商人要到n个不同的城市去推销商品，已知两两城市之间的距离，选择一条路径使商人经过所有城市并回到起点，使路程总长度最短。

搭建退火算法框架

• 我们先事先准备好实际问题所需要的数据，比如两两城市之间的距离矩阵等。
• 除此之外还要初始化退火算法的各种初始值，搭建温度循环和Markov链循环两个循环框架。对于解x和最优结果E，要设置三个参数new、current、best，分别用来存储当前值，上一轮的结果，历史最优解。

a = 0.99; T0 = 97; Tf = 3; T = T0; Markov_length = 10000;
coordinates_x = [565 25 345 945 845 880 25 525 580 650 1605 1220 1465 ...
    1530 845 725 145 415 510 560 300 520 480 835 975 1215 1320 1250 660 ...
    410 420 575 1150 700 685 685 770 795 720 760 475 95 875 700 555 830 ...
    1170 830 605 595 1340 1740];
coordinates_y = [575 185 750 685 655 660 230 1000 1175 1130 620 580 200 5 ...
    680 370 665 635 875 365 465  585 415 625 580 245 315 400 180 250 555 ...
    665 1160 580 595 610 610 645 635 650 960 260 920 500 815 485 65 610 ...
    625 360 725 245];
amount = size(coordinates_x,2);
dist_matrix = zeros(amount);
for i = 1:amount
    for j = 1:amount
        dist_matrix(i,j) = sqrt((coordinates_x(i)-coordinates_x(j)).^2 +...
            (coordinates_y(i)-coordinates_y(j)).^2);
    end
end
sol_new = 1:amount; sol_current = sol_new; sol_best = sol_new;
E_new = inf; E_current = inf; E_best = inf;

while T >= Tf
    for r = 1:Markov_length
    
    //填入具体程序
    
    end
    T = T.*a;
end

实现具体算法

对实际问题，核心是如何扰动。本题开始时初始化城市序号数列sol，按照默认顺序初始化可以满足初始值的E很大，然后扰动可以设计为任意交换两个城市的次序。根据sol我们可以转换成实际中的路线。

    ind1 = 0; ind2 = 0;
    while(ind1 == ind2)
        ind1 = randi(amount);
        ind2 = randi(amount);
    end
    tmp1 = sol_new(ind1);
    sol_new(ind1) = sol_new(ind2);
    sol_new(ind2) = tmp1;
    E_new = 0;
    for i = 1:(amount - 1)
        E_new = E_new + dist_matrix(sol_new(i),sol_new(i+1));
    end
    E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));
    if E_new < E_current
        E_current = E_new;
        sol_current = sol_new;
        if E_new < E_best
            E_best = E_new;
            sol_best = sol_new;
        end
    else
        if rand <exp(-(E_new - E_current)./T)
            E_current = E_new;
            sol_current = sol_new;
        else
            sol_new = sol_current;
        end
    end

完整程序

a = 0.99; T0 = 97; Tf = 3; T = T0; Markov_length = 10000;
coordinates_x = [565 25 345 945 845 880 25 525 580 650 1605 1220 1465 ...
    1530 845 725 145 415 510 560 300 520 480 835 975 1215 1320 1250 660 ...
    410 420 575 1150 700 685 685 770 795 720 760 475 95 875 700 555 830 ...
    1170 830 605 595 1340 1740];
coordinates_y = [575 185 750 685 655 660 230 1000 1175 1130 620 580 200 5 ...
    680 370 665 635 875 365 465  585 415 625 580 245 315 400 180 250 555 ...
    665 1160 580 595 610 610 645 635 650 960 260 920 500 815 485 65 610 ...
    625 360 725 245];
amount = size(coordinates_x,2);
dist_matrix = zeros(amount);
for i = 1:amount
    for j = 1:amount
        dist_matrix(i,j) = sqrt((coordinates_x(i)-coordinates_x(j)).^2 +...
            (coordinates_y(i)-coordinates_y(j)).^2);
    end
end
sol_new = 1:amount; sol_current = sol_new; sol_best = sol_new;
E_new = inf; E_current = inf; E_best = inf;

while T >= Tf
    for r = 1:Markov_length
        ind1 = 0; ind2 = 0;
        while(ind1 == ind2)
            ind1 = randi(amount);
            ind2 = randi(amount);
        end
        tmp1 = sol_new(ind1);
        sol_new(ind1) = sol_new(ind2);
        sol_new(ind2) = tmp1;
        E_new = 0;
        for i = 1:(amount - 1)
            E_new = E_new + dist_matrix(sol_new(i),sol_new(i+1));
        end
        E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));
        if E_new < E_current
            E_current = E_new;
            sol_current = sol_new;
            if E_new < E_best
                E_best = E_new;
                sol_best = sol_new;
            end
        else
            if rand <exp(-(E_new - E_current)./T)
                E_current = E_new;
                sol_current = sol_new;
            else
                sol_new = sol_current;
            end
        end
    end
    T = T.*a;
end

讨论退火参数对结果的影响

将上面的程序写成函数，将四个重要参数定义为函数的参数，即初始温度、终止温度、衰减因子、Markov链长度。
用控制变量法，分别对其中一个参数进行变动，其余参数保持原程序的经典值。为避免偶然性，每种情况计算十次取平均值。结果如下图所示：




从图中我们可以看出，参数在一定范围内有正常的波动性，但是如果偏离理想值太远，会严重影响结果。