模拟退火算法简介，以解决旅行商问题为例（基于Matlab)

模拟退火算法

基本概念

借由材料统计力学中的研究成果，通过模拟粒子在降温过程中的每个温度下达到热平衡的模式，最终得到处于低能状态的晶体（较为可靠的“全局最有解”）。

算法

基本参量

参数名称	含义
$T_0$	初始状态温度
$T$	当前状态温度
$T_{end}$	退火末态温度
$k_{max}$	平衡过程最大迭代次数
$\alpha$	$降温率T_{i+1}=\alpha T_i$
$x_0$	初始值
$x_i$	第i次迭代后状态
$x_{new}$	待定状态
$f(x)$	目标函数值（类比能量）
$neighbor(x)$	获得当前状态领域内的某个状态值
$P_{x_i\to x_{i+1}}$	当前状态转移至下一状态的概率

算法伪代码

\qquad\qquad
$T=T_0$
$x=x_0$
$while$ $T>T_{end}$
$fork=1to{k}_{m}ax$
$x_{new}=neighbor(x_i)$
$iff(x_{new})<f(x_i)$
$x_{i+1}=x_{new}$
$else$
$P_{x_i\to x_{i+1}}=exp(\frac{-[f(x_{new})-f(x_i)]}{T})$
$p=P_{x_i\to x_{i+1}}$
$ifp>r\leftarrow r\in U[0,1]$
$x_{i+1}=x_{new}$
$else$
$x_{i+1}=x_i$
$end$
$end$
$end$
$T=\alpha T$
$end$

注意

模拟退火算法的核心在于其以一定的概率去接受一个非理想状态下的解，借此以避免深陷于某个局部最小值。
整个优化过程是一个不断寻找新解和缓慢降温的交替过程，若温度下降十分缓慢，且在每个温度都有足够次数的温度转移，使之在每一个温度下达到热平衡，则全局最优解将一定会被找到。

1. 我们须综合考虑降温速度与解的质量的关系，在两者之间采取一种折中。若降温系数较小，降温过程迅速，则会有较大概率无法找到全局最优解。而若降温系数较大，降温过程过于缓慢，则其求解速度相较于其他搜索算法将不有速度优势。
2. 我们要确定在每一温度下状态转换的结束准则。如考虑当连续若干次的转换过程没有使状态发生变化是结束该温度下的状态转换。
3. 选择初始温度和确定某个可行解的邻域的方法也要恰当。

应用举例 —— 旅行商问题

旅行商问题(Traveling Salesman Problem, TSP)是一个经典的组合优化问题。经典的TSP可以描述为：一个商品推销员要去若干个城市推销商品，该推销员从一个城市出发，需要经过所有城市后，回到出发地。应如何选择行进路线，以使总的行程最短。

模型建立

内涵	模型
n个城市	n个节点，存入一个n维数组
访问城市顺序	n维数组中节点出现的先后顺序
邻域内新解的产生	任选序号i,j，调换这两者的顺序
目标函数	依次连接各节点产生的折线的长度

Matlab代码求解

数据准备

随机生成n个城市的坐标，并存入xls文档供后续使用
这里我随机生成了30个城市，范围在$[0,5000]\times [0,5000]$

%随机生成n个城市位置
%x,y range [0,5000]
clear,clc;
n = 30;
pos = randi([0,5000],n,2);
xs = pos(:,1);
ys = pos(:,2);
scatter(xs,ys,10,'filled');
%将数据写入文件
[stat,mes]=xlswrite('pos_data',pos);

先用蒙特卡洛法求得一个较好的初始解

function x_0 = get_initial_val(pos) %用蒙特卡洛法获得一个较好的初始值
k = 500;
for i = 1:k 
    new = shuffle(pos);
    if get_tour_length(new)<get_tour_length(pos)
        pos = new;
    end
end
x_0 = pos;
end

function out = shuffle(pos)
cur = 1;
len = length(pos);
while (cur < len/2)
    i = cur;
    j = len - cur + 1;
    change = randi([i,j]);
    pos = swch(pos,i,change);
    cur = cur + 1;
end
out = pos;
end

function res = swch(pos,a,b)
res = pos;
len = length(pos);
if len < a || len < b 
    return
end
if a < 1 || b < 1 || a==b
    return
end
temp = pos(a,:);
pos(a,:) = pos(b,:);
pos(b,:) = temp;
res = pos;
end

模拟退火算法主体

注：可通过调整初始温度，最终温度，温度衰减系数，最大迭代次数来改变结果的精度以及程序的运行速度

%SA算法解决旅行商问题
clear,clc;
%读取城市位置
pos = readmatrix('pos_data.xls');

T = 5000.0; %初始温度
T_end = 1e-8; %最终温度
alpha = 0.99; %温度衰减系数
k_max = 1000; %同一温度下迭代次数

pos = get_initial_val(pos); %通过蒙特卡洛法获得恰当的初始值
x0 = pos;
tic;
while T > T_end
    for k = 1:k_max
        new = get_neighbor(pos);
        if f(new)<=f(pos)
            pos = new;
        else
            prob = exp((f(pos)-f(new))/T);
            judge = rand;
            if prob > judge
                pos = new;
            end
        end
    end
    T = alpha*T;
    if T < 1000 && T> 990
        x1 = pos;
    elseif T< 100 && T> 99
        x2 = pos;
    end
        
end

display(f(pos));
toc;

subplot(2,2,1);
plot_tour(x0)
subplot(2,2,2);
plot_tour(x1);
subplot(2,2,3);
plot_tour(x2);
subplot(2,2,4);
plot_tour(pos);

function len = get_tour_length(pos) %计算路径总长
len = 0;
n = length(pos);
if n <= 1
    return
end
for i = 1:n-1
    x1 = pos(i,:);
    j = i+1;
    x2 = pos(j,:);
    len = len + get_distance(x1,x2);
end
len = len + get_distance(pos(1,:),pos(end,:));
end

function distance = get_distance(x1,x2)
dx=x1(1) - x2(1);
dy=x1(2) - x2(2);
distance = sqrt(dx^2+dy^2);
end

function out = f(pos)
out = get_tour_length(pos);
end

function plot_tour(pos)
xs = [pos(:,1);pos(1,1)];
ys = [pos(:,2);pos(1,2)];
plot(xs,ys,'-*');
end

function out = get_neighbor(pos)
i=1;
j=1;
while i == j
ind = randi([1,length(pos)],1,2);
i = ind(1);
j = ind(2);
end
temp = pos(i,:);
pos(i,:) = pos(j,:);
pos(j,:) = temp;
out = pos;
end

最终结果

旅行商周游一圈最短距离约为 $2.8659\times 10^4$
程序用时 36.355s

可视化

左上：由蒙特卡洛算法得出的初始解
右上：当T在1000附近时的解路线
左下：当T在100附近时的解路线
右下：最终解