[MCM] 2017研究生数学建模竞赛A题 3架飞机 TSP 求总路径最小

%% 数据准备
% 清空环境变量
clear all;
clc;

citys = xlsread('c:\Users\clemente\Desktop\terminal.xlsx','B2:C73');
ab1= [29 29 45.5 92 92 29]; % 五边形各点横坐标 最后回到第一个点 形成闭环
cd1= [125 58 57 85 125 125]; % 五边形各点纵坐标 最后回到第一个点 形成闭环
ab2= [45.5 92 92 45.5]; % 三角形各点横坐标 最后回到第一个点 形成闭环
cd2= [57 85 29 57]; % 三角形各点纵坐标 最后回到第一个点 形成闭环
ab3= [29 45.5 92 92 29 29]; % 五边形各点横坐标 最后回到第一个点 形成闭环
cd3= [58 57 29 12.5 13 58]; % 五边形各点纵坐标 最后回到第一个点 形成闭环

x=citys(:,1);
y=citys(:,2);
in1=inpolygon(x,y,ab1,cd1);
in2=inpolygon(x,y,ab2,cd2);
in3=inpolygon(x,y,ab3,cd3);

figure(1);
subplot(1,3,1);
plot(ab1,cd1,x(in1),y(in1),'ro',x(~in1),y(~in1),'b<')
subplot(1,3,2);
plot(ab2,cd2,x(in2),y(in2),'ro',x(~in2),y(~in2),'b<')
subplot(1,3,3);
plot(ab3,cd3,x(in3),y(in3),'ro',x(~in3),y(~in3),'b<')

city1=citys(in1,:);
city2=citys(in2,:);
city3=citys(in3,:);
basepoint=[110 0]; %无人机出发 基地的坐标
city1=[basepoint;city1;basepoint];  %基地加入到第一个聚类区域
city2=[basepoint;city2;basepoint];  %基地加入到第二个聚类区域
city3=[basepoint;city3;basepoint];  %基地加入到第三个聚类区域
%% 第一个聚类群

% 程序运行计时开始
t0 = clock;
%% 计算城市间距离
n = size(city1,1) -1; %城市数 起点,终点同一个要减1
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( city1(i,:) - city1(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = 27+1;                     % 蚂蚁数量最好接近城市数量的1.5倍 +1是因为新加入了基地
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = 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);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵;
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'']);
 
%% 绘图
figure(2)
plot([ city1(Shortest_Route,1);city1(Shortest_Route(1),1) ], [ city1(Shortest_Route,2);city1(Shortest_Route(1),2) ], 'k.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(city1,1)
    text(city1(i,1),city1(i,2),['  ' num2str(i)]);
end
text(city1(1,1),city1(1,2),'        起飞基地');

xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);

figure(3); %收敛图
subplot(1,3,1);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观


%% 第二个聚类群

% 程序运行计时开始
t0 = clock;
%% 计算城市间距离
n = size(city2,1)-1; %城市数 起点,终点同一个要减1
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( city2(i,:) - city2(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = 17+1;                     % 蚂蚁数量最好接近城市数量的1.5倍
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = 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);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵;
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'']);
 
%% 绘图
figure(2)
hold on;
plot([ city2(Shortest_Route,1);city2(Shortest_Route(1),1) ], [ city2(Shortest_Route,2);city2(Shortest_Route(1),2) ], 'r.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(city2,1)
    text(city2(i,1),city2(i,2),['  ' num2str(i)]);
end

xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);
hold off

figure(3); %收敛图
subplot(1,3,2);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观

%% 第三个聚类群

% 程序运行计时开始
t0 = clock;
%% 计算城市间距离
n = size(city3,1)-1; %城市数 起点,终点同一个要减1
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( city3(i,:) - city3(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = 27+1;                     % 蚂蚁数量最好接近城市数量的1.5倍 
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = 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);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵;
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'']);
 
%% 绘图
figure(2)
hold on;
plot([ city3(Shortest_Route,1);city3(Shortest_Route(1),1) ], [ city3(Shortest_Route,2);city3(Shortest_Route(1),2) ], 'b.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(city3,1)
    text(city3(i,1),city3(i,2),['  ' num2str(i)]);
end


xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);
hold off
axis([0 115 0 140]); %使得figure(2) 画布尺寸为100*120

figure(3); %收敛图
subplot(1,3,3);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观
多人 mTSP

代码改进

三处m值 改为

m = size(city1,1) -1;

m = size(city2,1) -1;

m = size(city3,1) -1;

避免维度错误

首先通过kmean(python)进行聚类,按无人机数量要求进行分类

数据附件下载

将聚类好的位置数据分别存入city1,city2,city3中分别应用蚁群算法进行TSP运算

 

缺点是由于没有考虑起点的距离 这里聚类应用不是特别好 导致上部分的路程较大 几乎是蓝红的两倍

最短距离:412.8601

最短路径:4 6 23 3 2 27 24 25 10 9 15 8 7 5 12 13 21 22 26 19 14 20 11 18 17 28 16 1 4

 

 

最短距离:209.7699

最短路径:7 18 15 17 11 16 12 14 19 2 3 5 4 8 9 6 10 13 1 7

 

 

最短距离:226.1174

最短路径:10 14 21 13 7 22 28 9 16 12 27 2 3 25 6 4 18 15 5 24 8 19 11 20 26 17 23 1 10

 

需要改进

 

%% 数据准备
% 清空环境变量
clear all;
clc;

citys = xlsread('c:\Users\clemente\Desktop\terminal.xlsx','B2:C73');
ab1= [91 27.6 27.08 91]; % 五边形各点横坐标 最后回到第一个点 形成闭环
cd1= [13  10.1 107.77 13]; % 五边形各点纵坐标 最后回到第一个点 形成闭环
ab2= [91 27.08 53.7  91]; % 三角形各点横坐标 最后回到第一个点 形成闭环
cd2= [13 107.77 129.4  13]; % 三角形各点纵坐标 最后回到第一个点 形成闭环
ab3= [91 53.7  93.8 91]; % 五边形各点横坐标 最后回到第一个点 形成闭环
cd3= [13 129.4 129.1 13]; % 五边形各点纵坐标 最后回到第一个点 形成闭环

x=citys(:,1);
y=citys(:,2);
in1=inpolygon(x,y,ab1,cd1);
in2=inpolygon(x,y,ab2,cd2);
in3=inpolygon(x,y,ab3,cd3);

figure(1);
subplot(1,3,1);
plot(ab1,cd1,x(in1),y(in1),'ro',x(~in1),y(~in1),'b<')
subplot(1,3,2);
plot(ab2,cd2,x(in2),y(in2),'ro',x(~in2),y(~in2),'b<')
subplot(1,3,3);
plot(ab3,cd3,x(in3),y(in3),'ro',x(~in3),y(~in3),'b<')

city1=citys(in1,:);
city2=citys(in2,:);
city3=citys(in3,:);
basepoint=[110 0]; %无人机出发 基地的坐标
city1=[basepoint;city1;basepoint];  %基地加入到第一个聚类区域
city2=[basepoint;city2;basepoint];  %基地加入到第二个聚类区域
city3=[basepoint;city3;basepoint];  %基地加入到第三个聚类区域
%% 第一个聚类群

% 程序运行计时开始
t0 = clock;
%% 计算城市间距离
n = size(city1,1) -1; %城市数 起点,终点同一个要减1
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( city1(i,:) - city1(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = size(city1,1) -1;                     % 蚂蚁数量最好接近城市数量的1.5倍 +1是因为新加入了基地
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = 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);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵;
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'']);
 
%% 绘图
figure(2)
plot([ city1(Shortest_Route,1);city1(Shortest_Route(1),1) ], [ city1(Shortest_Route,2);city1(Shortest_Route(1),2) ], 'k.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(city1,1)
    text(city1(i,1),city1(i,2),['  ' num2str(i)]);
end
text(city1(1,1),city1(1,2),'        起飞基地');

xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);

figure(3); %收敛图
subplot(1,3,1);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观


%% 第二个聚类群

% 程序运行计时开始
t0 = clock;
%% 计算城市间距离
n = size(city2,1)-1; %城市数 起点,终点同一个要减1
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( city2(i,:) - city2(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = size(city2,1) -1;                     % 蚂蚁数量最好接近城市数量的1.5倍
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = 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);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵;
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'']);
 
%% 绘图
figure(2)
hold on;
plot([ city2(Shortest_Route,1);city2(Shortest_Route(1),1) ], [ city2(Shortest_Route,2);city2(Shortest_Route(1),2) ], 'r.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(city2,1)
    text(city2(i,1),city2(i,2),['  ' num2str(i)]);
end

xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);
hold off

figure(3); %收敛图
subplot(1,3,2);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观

%% 第三个聚类群

% 程序运行计时开始
t0 = clock;
%% 计算城市间距离
n = size(city3,1)-1; %城市数 起点,终点同一个要减1
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( city3(i,:) - city3(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = size(city3,1) -1;                     % 蚂蚁数量最好接近城市数量的1.5倍 
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = 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);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵;
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'']);
 
%% 绘图
figure(2)
hold on;
plot([ city3(Shortest_Route,1);city3(Shortest_Route(1),1) ], [ city3(Shortest_Route,2);city3(Shortest_Route(1),2) ], 'b.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(city3,1)
    text(city3(i,1),city3(i,2),['  ' num2str(i)]);
end


xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);
hold off
axis([0 115 0 140]); %使得figure(2) 画布尺寸为100*120

figure(3); %收敛图
subplot(1,3,3);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观
采用扇形分类法 改进了m

 

最短距离:350.3253
最短路径:15 25 18 13 26 17 32 5 7 33 30 6 8 11 4 2 3 9 22 19 10 28 16 14 23 24 31 12 29 20 34 21 27 1 15
收敛迭代次数:204
程序执行时间:5.667秒
最短距离:319.4375
最短路径:5 2 3 4 8 7 14 21 9 10 16 22 12 18 23 11 17 19 1 13 15 6 20 5
收敛迭代次数:2
程序执行时间:2.288秒
最短距离:301.3621
最短路径:17 16 14 5 4 15 3 7 8 6 10 11 9 12 18 13 1 2 17
收敛迭代次数:256
程序执行时间:1.324秒

 

还是与最优化差距较大

下一步考虑改进 近距离的点合并

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值