写在前面：由于数学建模需要，最近学了启发式算法，个人感觉不难，就是与概率统计挂钩。另外在测试的时候需要频繁调参数，这就要看运气了…。 本文中代码来自中国科学院力学研究所–周吕文。这篇博文的目的是在记录学习心得的基础上给各位有缘人提供帮助。

算法思路易懂，我主要讲述代码含义（PS:注释都是我自己写的，Ahhhh）

一.思路

模拟退火首先从某个初始候选解开始，当温度大于0时执行循环。

在循环中，通过随机扰动产生一个新的解，然后求得新解和原解之间的能量差，如果差小于0，则采用新解作为当前解。

如果差大于0，则采用一个当前温度与能量差成比例的概率来选择是否接受新解。温度越低，接受的概率越小，差值越大，同样接受概率越小。

是否接受的概率用此公式计算：p=exp(-ΔE/T)。这里ΔE为新解与原解的差，T为当前的温度。

由于温度随迭代次数逐渐降低，因此获得一个较差的解的概率较小。

典型的模拟退火算法还使用了蒙特卡洛循环，在温度降低之前，通过多次迭代来找到当前温度下比较好的解。

二.以TSP为例详解退火算法

问题：假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。

1. 构建城市距离矩阵
2. 生成初始解
3. 计算初始解的值（即目标函数值）
4. 在温度阈值的基础上进行循环迭代计算（按照蒙特卡洛方法，温度降低之前多次实验，即设置一个参数，当条件达到时，就会进行降温，之后继续按照蒙特卡洛方法实验），每次计算都要进行解的扰动（产生新解），即进行尝试计算（扰动的方式有多种），计算新解的值并于原解比较，如果差小于0，则采用新解作为当前解，如果差大于0，则采用一个当前温度与能量差成比例的概率来选择是否接受新解。

%模拟退火算法
clear;clc;

load china; % geographic information
plotcities(province, border, city); % 画地图

numberofcities = length(city);      % 城市数目
% distance matrix: dis(i,j) is the distance between city i and j.
dis = distancematrix(city);   %城市距离矩阵


temperature = 1000;                 % 初始温度
cooling_rate = 0.94;                % 冷却率
iterations = 1;                     % Initialize the iteration number. 迭代数目（等温步长）

% Initialize random number generator with "seed". 
rand('seed',0);                    

% 随机选出初始解
route = randperm(numberofcities);
% 初始解的总距离
previous_distance = totaldistance(route,dis);

% 这是用于在100次迭代后冷却当前温度的标志
temperature_iterations = 1;
% 这是一个标志，用于在200次迭代后绘制当前路线。
plot_iterations = 1;

% plot the current route
plotroute(city, route, previous_distance, temperature);

while 1.0 < temperature
    %随机生成相邻解决方案
    temp_route = perturb(route,'reverse');
    % 计算临时路线的总距离
    current_distance = totaldistance(temp_route, dis);
    %计算距离变化
    diff = current_distance - previous_distance;
    
    % Metropolis Algorithm
    if (diff < 0) || (rand < exp(-diff/(temperature)))
        route = temp_route;         %接受新解
        previous_distance = current_distance;
        
        % 更新迭代
        temperature_iterations = temperature_iterations + 1;
        plot_iterations = plot_iterations + 1;
        iterations = iterations + 1;
    end
    
    % 每迭代100次降温
    if temperature_iterations >= 100
       temperature = cooling_rate*temperature;
       temperature_iterations = 0;
    end
    
    % 每200次迭代绘制当前路线
    if plot_iterations >= 200
       plotroute(city, route, previous_distance,temperature);
       plot_iterations = 0;
    end
end

% 绘制最终解决方案
plotroute(city, route, previous_distance,temperature);

代码片段详解：

1. numberofcities = length(city);      % 城市数目
   % distance matrix: dis(i,j) is the distance between city i and j.
   dis = distancematrix(city);   %城市距离矩阵

   首先获得城市数目，在用自己编写的distancematrix函数计算城市距离矩阵（计算方式由题意确定）

2. temperature = 1000;                 % 初始温度
   cooling_rate = 0.94;                % 冷却率
   iterations = 1;                     % Initialize the iteration number. 
   temperature_iterations = 1;

   定义算法参数，初始温度为1000°C。冷却率为0.94，即T(新)=T(原)*0.94,大部分都是以这种方式模拟降温，当然你也可以换成自己想要的方式。temperature_iterations初始化为1就是按蒙特卡洛方法进行循环，也可理解为等温步长，即在相同温度下筛选解，当筛选到你预定的参数时就会降温，接着在循环往复。

3. % 随机选出初始解
   route = randperm(numberofcities);
   % 初始解的总距离
   previous_distance = totaldistance(route,dis);

   不管是退火算法还是遗传算法都是需要初始解的，只不过遗传算法是初始解集，即种群的概念（遗传算法见我的博文）。route就是产生的初始解，即旅行城市的序号排序。用totaldistance函数计算这个解的总距离。

4. while 1.0 < temperature

   当温度低于1.0°C时就停止计算，也意味着已经计算出（按当前参数）最优解。

5.  %随机生成相邻解决方案
       temp_route = perturb(route,'reverse');
       % 计算临时路线的总距离
       current_distance = totaldistance(temp_route, dis);
       %计算距离变化
       diff = current_distance - previous_distance;

   扰动初始解，产生邻解（邻解意味着每次的扰动不能太大，否则计算不太精确，容易跳过最优解。扰动的方式有多种！）。计算这个新解的总距离，并计算新解和原解的差值diff

6.  if (diff < 0) || (rand < exp(-diff/(temperature)))
           route = temp_route;         %接受新解
           previous_distance = current_distance;
           
           % 更新迭代
           temperature_iterations = temperature_iterations + 1;
           plot_iterations = plot_iterations + 1;
           iterations = iterations + 1;
       end

   如果diff<0（即新解的值比原解小）或满足Metropolis准则——以概率接受新状态，则接受这个解。并更新当前的总距离和把原解替换为这个新解。temperature_iterations自增一(这里是每接受一个新解就会增1，当然你也可以按照自己的意愿设置它的自增规则，但它一定是要有的)

7.   % 每迭代100次降温
       if temperature_iterations >= 100
          temperature = cooling_rate*temperature;
          temperature_iterations = 0;
       end

   当相同温度下接受新解的个数累积到100次时进行降温，并把temperature_iterations置0.

8. 温度达到阈值，程序终止，得到最优解。

扰动函数：不解释，就是2种扰动方式

function route = perturb(route_old, method)
% PERTURB
% route = PERTURB(route_old, method) generate randomly a neighbouring route by
% perturb old route. perturb methods:
%                        ___________            ___________         
%     1. reverse:   [1 2 3 4 5 6 7 8 9] -> [1 2 8 7 6 5 4 3 9]
%                        _         _            _         _
%     2. swap:      [1 2 3 4 5 6 7 8 9] -> [1 2 8 4 5 6 7 3 9]

route = route_old;
numbercities = length(route);
city1 = ceil(numbercities*rand);
city2 = ceil(numbercities*rand);
switch method
    case 'reverse'
        citymin = min(city1,city2);
        citymax = max(city1,city2);
        route(citymin:citymax) = route(citymax:-1:citymin);
    case 'swap'
        route([city1, city2]) = route([city2, city1]);
end

计算总距离函数：

function d = totaldistance(route, dis)
% TOTALDISTANCE
% d = TOTALDISTANCE(route, dis) calculate total distance of a route with
% the distance matrix dis.

d = dis(route(end),route(1)); % closed path
for k = 1:length(route)-1
    i = route(k);
    j = route(k+1);
    d = d + dis(i,j);
end

以上就是算法的核心，其他的代码就是绘制图形。

运行结果如下：

完整代码：

https://github.com/Galgaddott/-

PS:最优解是当前参数下的最优解，可能是局部最优解，因为算法都是有局限的，所以需要自己多次调整参数，直到找到另你满意的解。