遗传，模拟退火，蚁群三个算法求解TSP的对比 - 图文(4)

模拟退火:

clear clc

a = 0.99; % 温度衰减函数的参数 t0 = 97; tf = 3; t = t0; Markov_length = 10000; % Markov链长度

citys=[1 1304 2312;2 3639 1315;3 4177 2244;4 3721 1399;5 3488 1535;6 3326 1556;7 3238 1229;8 4196 1004;9 4312 790;10 4386 570;11 3007 1970;12 2562 1756;13 2788 1491;14 2381 1676;15 1332 695;16 3715 1678;17 3918 2179;18 4061 2370;19 3780 2212;20 3676 2578;21 4029 2838;22 4263 2931;23 3429 1908;24 3507 2367;25 3394 2643;26 3439 3201;27 2935 3240;28 3140 3550;29 2545 2357;30 2778 2826;31 2370 2975]; citys(:,1) = [];

amount = size(citys,1); % 城市的数目 % 通过向量化的方法计算距离矩阵 dist_matrix = zeros(amount, amount);

coor_x_tmp1 =citys(:,1) * ones(1,amount); coor_x_tmp2 = coor_x_tmp1';

coor_y_tmp1 = citys(:,2) * ones(1,amount); coor_y_tmp2 = coor_y_tmp1';

dist_matrix = sqrt((coor_x_tmp1-coor_x_tmp2).^2 + ... (coor_y_tmp1-coor_y_tmp2).^2);

sol_new = 1:amount; % 产生初始解

% sol_new是每次产生的新解；sol_current是当前解；sol_best是冷却中的最好解； E_current = inf; % E_current是当前解对应的回路距离； E_best = inf; % E_best是最优解的 % E_new是新解的回路距离；

sol_current = sol_new; sol_best = sol_new; p = 1;

while t>=tf

for r=1:Markov_length % Markov链长度 % 产生随机扰动

if (rand < 0.5) % 随机决定是进行两交换还是三交换 % 两交换

ind1 = 0; ind2 = 0; while (ind1 == ind2)

ind1 = ceil(rand.*amount); ind2 = ceil(rand.*amount); end

tmp1 = sol_new(ind1);

sol_new(ind1) = sol_new(ind2); sol_new(ind2) = tmp1; else

% 三交换

ind1 = 0; ind2 = 0; ind3 = 0;

while (ind1 == ind2) || (ind1 == ind3) ... || (ind2 == ind3) || (abs(ind1-ind2) == 1) ind1 = ceil(rand.*amount); ind2 = ceil(rand.*amount); ind3 = ceil(rand.*amount); end

tmp1 = ind1;tmp2 = ind2;tmp3 = ind3; % 确保ind1 < ind2 < ind3

if (ind1 < ind2) && (ind2 < ind3) ;

elseif (ind1 < ind3) && (ind3 < ind2) ind2 = tmp3;ind3 = tmp2;

elseif (ind2 < ind1) && (ind1 < ind3) ind1 = tmp2;ind2 = tmp1;

elseif (ind2 < ind3) && (ind3 < ind1)

ind1 = tmp2;ind2 = tmp3; ind3 = tmp1; elseif (ind3 < ind1) && (ind1 < ind2) ind1 = tmp3;ind2 = tmp1; ind3 = tmp2; elseif (ind3 < ind2) && (ind2 < ind1) ind1 = tmp3;ind2 = tmp2; ind3 = tmp1; end

tmplist1 = sol_new((ind1+1):(ind2-1)); sol_new((ind1+1):(ind1+ind3-ind2+1)) = ... sol_new((ind2):(ind3));

sol_new((ind1+ind3-ind2+2):ind3) = ... tmplist1; end

%检查是否满足约束

% 计算目标函数值（即内能） 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; % 控制参数t（温度）减少为原来的a倍 end

disp('最优解为：') disp(sol_best)

disp('最短距离：') disp(E_best)

程序运行结果为:

蚁群算法:

clear all clc

%% 导入数据

citys=[1304 2312;3639 1315;4177 2244;3721 1399;3488 1535;3326 1556;3238 1229;4196 1004;4312 790;4386 570;3007 1970;2562 1756;2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2367;3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975]; %% 计算城市间相互距离 n = size(citys,1); D = zeros(n,n); for i = 1:n for j = 1:n if i ~= j

D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2)); else

D(i,j) = 1e-4; end end end

%% 初始化参数

m = 50; % 蚂蚁数量

alpha = 1; % 信息素重要程度因子 beta = 5; % 启发函数重要程度因子 rho = 0.1; % 信息素挥发因子 Q = 1; % 常系数 Eta = 1./D; % 启发函数 Tau = ones(n,n); % 信息素矩阵 Table = zeros(m,n); % 路径记录表 iter = 1; % 迭代次数初值 iter_max = 200; % 最大迭代次数 Route_best = zeros(iter_max,n); % 各代最佳路径 Length_best = zeros(iter_max,1); % 各代最佳路径的长度 Length_ave = zeros(iter_max,1); % 各代路径的平均长度

%% 迭代寻找最佳路径 while iter <= iter_max

% 随机产生各个蚂蚁的起点城市 start = zeros(m,1); for i = 1:m

temp = randperm(n); start(i) = temp(1); end

Table(:,1) = start; % 构建解空间 citys_index = 1:n;

% 逐个蚂蚁路径选择 for i = 1:m

% 逐个城市路径选择 for j = 2:n

tabu = Table(i,1:(j - 1)); % 已访问的城市集合(禁忌表) allow_index = ~ismember(citys_index,tabu);

allow = citys_index(allow_index); % 待访问的城市集合 P = allow;

% 计算城市间转移概率 for k = 1:length(allow)

P(k) = Tau(tabu(end),allow(k))^alpha * Eta(tabu(end),allow(k))^beta; end

P = P/sum(P);

% 轮盘赌法选择下一个访问城市 Pc = cumsum(P);

target_index = find(Pc >= rand); target = allow(target_index(1));