快速了解模拟退火算法-基于matlab实现

 

 

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

算法思路易懂,我主要讲述代码含义(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:最优解是当前参数下的最优解,可能是局部最优解,因为算法都是有局限的,所以需要自己多次调整参数,直到找到另你满意的解。

 

展开阅读全文
©️2020 CSDN 皮肤主题: 大白 设计师: CSDN官方博客 返回首页
实付0元
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值