码迷,mamicode.com
首页 > 编程语言 > 详细

模拟退火算法-旅行商问题-matlab实现

时间:2019-01-20 18:00:57      阅读:272      评论:0      收藏:0      [点我收藏+]

标签:load   stat   inf   count   mat   不同   end   之间   模拟退火   

整理一下数学建模会用到的算法,供比赛时候参考食用。

——————————————————————————————————————————

旅行商问题(TSP):

给定一系列城市和每对城市之间的距离,求解访问每一座城市一次并回到起始城市的最短回路。

它是组合优化中的一个NP困难问题,在运筹学和理论计算机科学中非常重要。

以下两个程序,在不同数据集合情况下表现有所差别,理论上第一个程序的解更为优化。

  1 clear
  2 clc
  3 a = 0.99;       %温度衰减函数的参数
  4 t0 = 97;        %初始温度
  5 tf = 3;         %终止温度
  6 t = t0;
  7 Markov_length = 10000;  %Markov链长度
  8  
  9 % load data.txt
 10 % x = data(:, 1:2:8); x = x(:);
 11 % y = data(:, 2:2:8); y = y(:);
 12 % data = [70,40;x, y];
 13 % coordinates = data;
 14 coordinates = [
 15     1    565.0   575.0; 2     25.0   185.0; 3    345.0   750.0;
 16     4    945.0   685.0; 5    845.0   655.0; 6    880.0   660.0;
 17     7     25.0   230.0; 8    525.0  1000.0; 9    580.0  1175.0;
 18     10   650.0  1130.0; 11  1605.0   620.0; 12  1220.0   580.0;
 19     13  1465.0   200.0; 14  1530.0     5.0; 15   845.0   680.0;
 20     16   725.0   370.0; 17   145.0   665.0; 18   415.0   635.0;
 21     19   510.0   875.0; 20   560.0   365.0; 21   300.0   465.0;
 22     22   520.0   585.0; 23   480.0   415.0; 24   835.0   625.0;
 23     25   975.0   580.0; 26  1215.0   245.0; 27  1320.0   315.0;
 24     28  1250.0   400.0; 29   660.0   180.0; 30   410.0   250.0;
 25     31   420.0   555.0; 32   575.0   665.0; 33  1150.0  1160.0;
 26     34   700.0   580.0; 35   685.0   595.0; 36   685.0   610.0;
 27     37   770.0   610.0; 38   795.0   645.0; 39   720.0   635.0;
 28     40   760.0   650.0; 41   475.0   960.0; 42    95.0   260.0;
 29     43   875.0   920.0; 44   700.0   500.0; 45   555.0   815.0;
 30     46   830.0   485.0; 47  1170.0    65.0; 48   830.0   610.0;
 31     49   605.0   625.0; 50   595.0   360.0; 51  1340.0   725.0;
 32     52  1740.0   245.0;
 33     ];
 34 coordinates(:,1) = [];
 35 amount = size(coordinates,1);        %城市的数目
 36 %通过向量化的方法计算距离矩阵
 37 dist_matrix = zeros(amount,amount);
 38 coor_x_tmp1 = coordinates(:,1) * ones(1,amount);
 39 coor_x_tmp2 = coor_x_tmp1;
 40 coor_y_tmp1 = coordinates(:,2) * ones(1,amount);
 41 coor_y_tmp2 = coor_y_tmp1;
 42 dist_matrix = sqrt((coor_x_tmp1 - coor_x_tmp2).^2 + (coor_y_tmp1 - coor_y_tmp2).^2);
 43 
 44 
 45 sol_new = 1:amount;         %产生初始解,sol_new是每次产生的新解
 46 sol_current = sol_new;      %sol_current是当前解
 47 sol_best = sol_new;         %sol_best是冷却中的最好解
 48 E_current = inf;            %E_current是当前解对应的回路距离
 49 E_best = inf;               %E_best是最优解
 50 p = 1;
 51 
 52 
 53 rand(state, sum(clock));
 54 
 55 for j = 1:10000
 56     sol_current = [randperm(amount)];
 57     E_current = 0;
 58     for i=1:(amount-1)
 59         E_current = E_current+dist_matrix(sol_current(i), sol_current(i+1));
 60     end
 61     if E_current<E_best
 62         sol_best = sol_current;
 63         E_best = E_current;
 64     end
 65 end
 66 
 67 
 68 
 69 
 70 while t >= tf
 71    for r = 1:Markov_length      %Markov链长度
 72     %产生随机扰动
 73     if(rand < 0.5)
 74         %两交换
 75         ind1 = 0;
 76         ind2 = 0;
 77         while(ind1 == ind2)
 78            ind1 = ceil(rand * amount);
 79            ind2 = ceil(rand * amount);
 80         end
 81         tmp1 = sol_new(ind1);
 82         sol_new(ind1) = sol_new(ind2);
 83         sol_new(ind2) = tmp1;
 84     else
 85         %三交换
 86         ind1 = 0;
 87         ind2 = 0;
 88         ind3 = 0;
 89         while( (ind1 == ind2) || (ind1 == ind3) || (ind2 == ind3) || (abs(ind1 -ind2) == 1) )
 90             ind1 = ceil(rand * amount);
 91             ind2 = ceil(rand * amount);
 92             ind3 = ceil(rand * amount);
 93         end
 94         tmp1 = ind1;
 95         tmp2 = ind2;
 96         tmp3 = ind3;
 97         %确保 ind1 < ind2 < ind3
 98         if(ind1 < ind2) && (ind2 < ind3);
 99         elseif(ind1 < ind3) && (ind3 < ind2)
100             ind1 = tmp1; ind2 = tmp3; ind3 = tmp2;
101         elseif(ind2 < ind1) && (ind1 < ind3)
102             ind1 = tmp2; ind2 = tmp1; ind3 = tmp3;
103         elseif(ind2 < ind3) && (ind3 < ind1)
104             ind1 = tmp2; ind2 = tmp3; ind3 = tmp1;
105         elseif(ind3 < ind1) && (ind1 < ind2)
106             ind1 = tmp3; ind2 = tmp1; ind3 = tmp2;
107         elseif(ind3 < ind2) && (ind2 < ind1)
108             ind1 = tmp3; ind2 = tmp2; ind3 = tmp1;
109         end
110         
111         tmplist1 = sol_new((ind1 + 1):(ind2 - 1));
112         sol_new((ind1 + 1):(ind1 + (ind3 - ind2 + 1) )) = sol_new((ind2):(ind3));
113         sol_new((ind1 + (ind3 - ind2 + 1) + 1):(ind3)) = tmplist1;
114     end
115     
116     %检查是否满足约束
117     
118     %计算目标函数值(即内能)
119     E_new = 0;
120     for i = 1:(amount - 1)
121         E_new = E_new + dist_matrix(sol_new(i),sol_new(i + 1));
122     end
123     %再算上从最后一个城市到第一个城市的距离
124     E_new = E_new + dist_matrix(sol_new(amount),sol_new(1));
125     
126     if E_new < E_current
127         E_current = E_new;
128         sol_current = sol_new;
129         if E_new < E_best
130             E_best = E_new;
131             sol_best = sol_new;
132         end
133     else
134         %若新解的目标函数值大于当前解,
135         %则仅以一定概率接受新解
136         if rand < exp(-(E_new - E_current) / t)
137             E_current = E_new;
138             sol_current = sol_new;
139         else
140             sol_new = sol_current;
141         end
142         
143     end
144    end
145 
146    t = t * a;      %控制参数t(温度)减少为原来的a倍
147 end
148 
149 E_best = E_best+dist_matrix(sol_best(end), sol_best(1));
150 
151 disp(最优解为:);
152 disp(sol_best);
153 disp(最短距离:);
154 disp(E_best);
155 
156 data1 = zeros(length(sol_best),2 );
157 for i = 1:length(sol_best)
158     data1(i, :) = coordinates(sol_best(1,i), :);
159 end
160 
161 data1 = [data1; coordinates(sol_best(1,1),:)];
162 
163 
164 figure
165 plot(coordinates(:,1), coordinates(:,2), *k, data1(:,1), data1(:, 2), r); 
166 title( [ 近似最短路径如下,路程为 , num2str( E_best ) ] ) ;

另一种根据《数学建模算法与应用—司守奎》中算法改编:

clc;
clear;
close all;

coordinates = [
    2     25.0   185.0; 3    345.0   750.0;
    4    945.0   685.0; 5    845.0   655.0; 6    880.0   660.0;
    7     25.0   230.0; 8    525.0  1000.0; 9    580.0  1175.0;
    10   650.0  1130.0; 11  1605.0   620.0; 12  1220.0   580.0;
    13  1465.0   200.0; 14  1530.0     5.0; 15   845.0   680.0;
    16   725.0   370.0; 17   145.0   665.0; 18   415.0   635.0;
    19   510.0   875.0; 20   560.0   365.0; 21   300.0   465.0;
    22   520.0   585.0; 23   480.0   415.0; 24   835.0   625.0;
    25   975.0   580.0; 26  1215.0   245.0; 27  1320.0   315.0;
    28  1250.0   400.0; 29   660.0   180.0; 30   410.0   250.0;
    31   420.0   555.0; 32   575.0   665.0; 33  1150.0  1160.0;
    34   700.0   580.0; 35   685.0   595.0; 36   685.0   610.0;
    37   770.0   610.0; 38   795.0   645.0; 39   720.0   635.0;
    40   760.0   650.0; 41   475.0   960.0; 42    95.0   260.0;
    43   875.0   920.0; 44   700.0   500.0; 45   555.0   815.0;
    46   830.0   485.0; 47  1170.0    65.0; 48   830.0   610.0;
    49   605.0   625.0; 50   595.0   360.0; 51  1340.0   725.0;
    52  1740.0   245.0;
    ];
coordinates(:,1) = [];
data = coordinates;


% 读取数据
% load data.txt;

% x = data(:, 1:2:8); x = x(:);
% y = data(:, 2:2:8); y = y(:);
x = data(:, 1);
y = data(:, 2);
start = [565.0   575.0];
data = [start; data;start];

% data = [start; x, y;start];
% data = data*pi/180;


% 计算距离的邻接表
count = length(data(:, 1));
d = zeros(count);
for i = 1:count-1
    for j = i+1:count
%         temp = cos(data(i, 1)-data(j,1))*cos(data(i,2))*cos(data(j,2))...
%             +sin(data(i,2))*sin(data(j,2));
         d(i,j)=( sum( ( data( i , : ) - data( j , : ) ) .^ 2 ) ) ^ 0.5 ;
%         d(i,j) = 6370*acos(temp);
    end
end
d =d + d;  % 对称  i到j==j到i



S0=[];          % 存储初值
Sum=inf;     % 存储总距离

rand(state, sum(clock));

% 求一个较为优化的解,作为初值
for j = 1:10000
    S = [1 1+randperm(count-2), count];
    temp = 0;
    for i=1:count-1
        temp = temp+d(S(i), S(i+1));
    end
    if temp<Sum
        S0 = S;
        Sum = temp;
    end
end

e = 0.1^40;   % 终止温度
L = 2000000;     % 最大迭代次数
at = 0.999999;     % 降温系数
T = 2;             % 初温

% 退火过程
for k = 1:L
    % 产生新解
    c =1+floor((count-1)*rand(1,2));
    
    c = sort(c);
    c1 = c(1);  c2 = c(2);
    if c1==1
        c1 = c1+1;
    end
    if c2==1
        c2 = c2+1;
    end
    % 计算代价函数值
    df = d(S0(c1-1), S0(c2))+d(S0(c1), S0(c2+1))-...
        (d(S0(c1-1), S0(c1))+d(S0(c2), S0(c2+1)));
    % 接受准则
    if df<0
        S0 = [S0(1: c1-1), S0(c2:-1:c1), S0(c2+1:count)];
        Sum = Sum+df;
    elseif exp(-df/T) > rand(1)
        S0 = [S0(1: c1-1), S0(c2:-1:c1), S0(c2+1:count)];
        Sum = Sum+df;
    end
    T = T*at;
    if T<e
        break;
    end
end

data1 = zeros(2, count);
% y1 = [start; x, y; start];
for i =1:count
   data1(:, i) = data(S0(1,i), :);
end

figure
plot(x, y, o, data1(1, :), data1(2, :), r);
title( [ 近似最短路径如下,路程为 , num2str( Sum ) ] ) ;
disp(Sum);
S0

 

模拟退火算法-旅行商问题-matlab实现

标签:load   stat   inf   count   mat   不同   end   之间   模拟退火   

原文地址:https://www.cnblogs.com/Dawn-bin/p/10295557.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!