智能优化算法解决TSP问题

Aircraft GNC

已于 2022-04-13 18:27:18 修改

阅读量1.6k

点赞数 4

分类专栏：智能优化算法文章标签：算法

于 2020-04-08 00:02:37 首次发布

本文链接：https://blog.csdn.net/qq_43614222/article/details/105377313

版权

智能优化算法专栏收录该内容

1 篇文章 0 订阅

订阅专栏

遗传算法

%% 1.清空环境变量
clear;
clc;

%% 2.导入数据
load citys_data.mat;          %数据集的变量名为city

%% 3.计算城市间相互距离
n=size(city,1);
D=zeros(n,n);
for i=1:n
    for j=1:n
        if i~=j
            D(i,j)=sqrt(sum((city(i,:)-city(j,:)).^2));
        else
            D(i,j)=1e-4;
        end
    end
end

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

%% 5.迭代寻找最佳路径
h = waitbar(0,' Time Loop ');
tic;
while iter<=itermax
    %% 5.1随机产生各个蚂蚁的起点城市
    start=zeros(m,1);
    for i=1:m
        temp=randperm(n);
        start(i,:)=temp(1);
    end
    Table(:,1)=start;
    citys_index=1:n;

    %% 5.2产生解空间（路径表）
    %逐个蚂蚁路径选择
    for i=1:m
        %逐个城市路径选择
        for j=2:n
            tabu = Table(i,:); %已访问城市
            allow_index=~ismember(citys_index,tabu);
            allow=citys_index(allow_index);
            P=allow;

            %% 5.3计算城市间转移概率
            for k=1:length(allow) 
                P(k)=Tau(Table(i,j-1),allow(k))^alpha * Eta(Table(i,j-1),allow(k))^beta;
            end

            P=P/sum(P);

            %% 5.4轮盘赌法选择下一个访问城市
            Pc=cumsum(P);
            target_index = find(Pc>=rand);
            target=allow(target_index(1));
            Table(i,j)=target;
        end
    end

    %% 5.5计算各个蚂蚁的路径距离
    Length=zeros(m,1);
    for i=1:m
        for j=1:n-1
            Length(i,1)=Length(i,1)+D(Table(i,j),Table(i,j+1));
        end
        Length(i,1)=Length(i,1)+D(Table(i,end),Table(i,1));
    end

    %% 5.6计算最短路径及平均距离
    if iter==1
        [min_Length,min_index]=min(Length);
        Length_best(iter)=min_Length;
        Route_best(iter,:)=Table(min_index,:);
        Length_ave(iter)=mean(Length);
    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-1)>min_Length
            Route_best(iter,:)=Table(min_index,:);
        else
            Route_best(iter,:)=Route_best(iter-1,:);
        end
    end

    %% 5.7更新信息素
    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);
        end
        Delta_Tau(Table(i,end),Table(i,1))=Delta_Tau(Table(i,end),Table(i,1)) + Q/Length(i);
    end
    Tau=(1-rho)*Tau+Delta_Tau;

    % 5.8迭代次数加1，清空路径表
    iter=iter+1;
    Table=zeros(m,n);
    waitbar(iter / itermax, h, sprintf('%s%%', num2str(100 * iter / itermax)));

end
close(h);
toc;

%% 6.结果显示
best_route=Route_best(end,:);
best_length=Length_best(end,:);
OutputPath(best_route);                     %列出最短路径
disp(['最短距离: ' num2str(best_length)]);  %列出最短距离

%% 7.绘图
figure(1);
plot([city(best_route,1);city(best_route(1),1)],[city(best_route,2);city(best_route(1),2)],'o-');
for i=1:size(city,1)
    text(city(i,1),city(i,2),[' ' num2str(i)]);    % 标出城市序号
end
text(city(best_route(1),1),city(best_route(1),2),'    起点');  % 标出起点城市
text(city(best_route(end),1),city(best_route(end),2),'    终点'); % 标出终点城市
grid on;
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['蚁群算法优化路径(最短距离：' num2str(best_length) ')']);

figure(2);
plot(1:itermax,Length_ave,'r:',1:itermax,Length_best,'b-');
grid on;
legend('平均距离','最短距离');
xlabel('迭代次数');
ylabel('距离');
title('各代最短距离与平均距离对比');

function OutputPath(S)
%% 输出路径函数
% 输入：S 路径
S = [S,S(1)];
N = length(S);
p = num2str(S(1));  % p是字符串
for i = 2:N
    p = [p,'―>',num2str(S(i))];
end
disp(p);

粒子群算法

%% 智能优化算法系列之粒子群算法
%%  参数说明----------------------------------------------------------
% c1                         个体学习因子
% c2                         社会学习因子
% w                          惯性因子
% m                          粒子数量
% pop                        粒子位置
% v                          粒子速度
% gen                        迭代计数器
% genmax                     迭代次数
% fitness                    适应度函数值
% Pbest                      个体极值路径
% Pbest_fitness              个体极值
% Gbest                      群体极值路径
% Gbest_fitness              群体极值
% Length_ave                 各代路径的平均长度
% ws                         惯性因子最大值
% we                         惯性因子最小值

%% 1.清空环境变量
clear;
clc;
tic;

%% 2.导入数据
load citys_data.mat;          %数据集的变量名为city

%% 3.计算城市间相互距离
n = size(city,1);
d = zeros(n,n);

for i=1:n
    for j=i+1:n
        d(i,j) = sqrt(sum((city(i,:) - city(j,:)).^2));
        d(j,i) = d(i,j);
    end
end

%% 4.初始化参数
c1 = 0.1;                         %个体学习因子
c2 = 0.075;                       %社会学习因子
w = 1;                            %惯性因子
m = 500;                          %粒子数量
pop = zeros(m,n);                 %粒子位置
v = zeros(m,n);                   %粒子速度
gen = 1;                          %迭代计数器
genmax = 2000;                    %迭代次数
fitness = zeros(m,1);             %适应度函数值
Pbest = zeros(m,n);               %个体极值路径
Pbest_fitness = zeros(m,1);       %个体极值
Gbest = zeros(genmax,n);          %群体极值路径
Gbest_fitness = zeros(genmax,1);  %群体极值
Length_ave = zeros(genmax,1);     %各代路径的平均长度
ws = 1;                           %惯性因子最大值
we = 0.8;                         %惯性因子最小值

%% 5.产生初始粒子
% 5.1随机产生粒子初始位置和速度
for i=1:m
    pop(i,:) = randperm(n);
    v(i,:) = randperm(n);
end

% 5.2计算粒子适应度函数值
for i=1:m
    for j=1:n-1
        fitness(i) = fitness(i) + d(pop(i,j),pop(i,j+1));
    end
    fitness(i) = fitness(i) + d(pop(i,end),pop(i,1));
end

% 5.3计算个体极值和群体极值
Pbest_fitness = fitness;
Pbest = pop;
[Gbest_fitness(1),min_index] = min(fitness);
Gbest(1,:) = pop(min_index,:);
Length_ave(1) = mean(fitness);

%% 6.迭代寻优
h = waitbar(0,' Time Loop '); %进度条
while gen<genmax
    %% 6.1更新迭代次数与惯性因子
    gen = gen + 1;
    w = ws - (ws-we) * (gen/genmax)^2;

    %% 6.2更新速度
    %个体极值修正部分
    change1 = position_minus_position(Pbest,pop);
    change1 = constant_times_velocity(c1,change1);
    %群体极值修正部分
    change2 = position_minus_position(repmat(Gbest(gen-1,:),m,1),pop);
    change2 = constant_times_velocity(c2,change2);
    %原速度部分
    v = constant_times_velocity(w,v);
    %修正速度
    for i=1:m
        for j=1:n
            if change1(i,j)~=0
                v(i,j) = change1(i,j);
            end
            if change2(i,j)~=0
                v(i,j) = change2(i,j);
            end
        end
    end

    %% 6.3更新位置
    pop = position_plus_velocity(pop,v);

    %% 6.4适应度函数值更新
    fitness = zeros(m,1);
    for i=1:m
        for j=1:n-1
            fitness(i) = fitness(i) + d(pop(i,j),pop(i,j+1));
        end
        fitness(i) = fitness(i) + d(pop(i,end),pop(i,1));
    end

    %% 6.5个体极值与群体极值更新
    for i=1:m
        if fitness(i)<Pbest_fitness(i)
            Pbest_fitness(i) = fitness(i);
            Pbest(i,:) = pop(i,:);
        end
    end

    [minvalue,min_index] = min(fitness);
    if minvalue<Gbest_fitness(gen-1)
        Gbest_fitness(gen) = minvalue;
        Gbest(gen,:) = pop(min_index,:);
    else
        Gbest_fitness(gen) = Gbest_fitness(gen-1);
        Gbest(gen,:) = Gbest(gen-1,:);
    end

    Length_ave(gen) = mean(fitness);
    waitbar(gen / genmax, h, sprintf('%s%%', num2str(100 * gen / genmax)));
end
close(h);
%% 7.结果显示
[Shortest_Length,Shortest_index] = min(Gbest_fitness);
Shortest_Route = Gbest(Shortest_index,:);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str([Shortest_Route Shortest_Route(1)])]);

%% 8.绘图
figure(1);
plot([city(Shortest_Route,1);city(Shortest_Route(1),1)],[city(Shortest_Route,2);city(Shortest_Route(1),2)],'o-');
grid on;
for i = 1:size(city,1)
    text(city(i,1),city(i,2),['   ' num2str(i)]);
end
text(city(Shortest_Route(1),1),city(Shortest_Route(1),2),'       起点');
text(city(Shortest_Route(end),1),city(Shortest_Route(end),2),'       终点');
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['粒子群算法优化路径(最短距离:' num2str(Shortest_Length) ')']);

figure(2);
plot(1:genmax,Gbest_fitness,'b',1:genmax,Length_ave,'r:');
grid on;
legend('最短距离','平均距离');
xlabel('迭代次数');
ylabel('距离');
title('各代最短距离与平均距离对比');
toc;

function pop = position_plus_velocity(pop,v)
%利用速度记录的交换序列进行位置修正
%6.3更新位置
for i=1:size(pop,1)
    for j=1:size(pop,2)
        if v(i,j)~=0
            temp = pop(i,j);
            pop(i,j) = pop(i,v(i,j));
            pop(i,v(i,j)) = temp;
        end
    end
end
end

function change = position_minus_position(best,pop)
%记录将pop变成best的交换序列
% 6.2更新速度
%个体极值修正部分
for i=1:size(best,1)
    for j=1:size(best,2)
        change(i,j) = find(pop(i,:)==best(i,j)); % return ture or false 1\0
        temp=pop(i,j);
        pop(i,j)=pop(i,change(i,j));
        pop(i,change(i,j))=temp;
    end
end
end

function change = constant_times_velocity(constant,change)
% 以一定概率保留交换序列
for i=1:size(change,1)
    for j=1:size(change,2)
        if (rand > constant) 
            change(i,j) = 0;
        end
    end
end
end

模拟退火算法

%% I. 清空环境变量
clear;
clc;
tic;
%% 2.导入数据
load citys_data.mat;          %数据集的变量名为city
 
%% III. 计算距离矩阵
D = Distance(city);  %计算距离矩阵
N = size(city,1);    %城市的个数
 
%% IV. 初始化参数
T0 = 10^(20);   % 初始温度，10的10次方！需要设定一个很大的温度。
T_end = 1e-50;   % 终止温度
L = 2;          % 各温度下的迭代次数
q = 0.9;        % 降温速率
Time = 132;
count = 0;                   %迭代计数
Obj = zeros(Time,1);         %目标值矩阵初始化
track = zeros(Time,N);       %每代的最优路线矩阵初始化
 
%% V. 随机产生一个初始路线
S1 = randperm(N);             % 返回一个大小为 N 的值为1->N的向量

%% VI. 迭代优化
while (T0 > T_end)
    %% 更新迭代次数
    count = count + 1;     
    temp = zeros(L,N+1);
    
    %% 1. 产生新解
    S2 = NewAnswer(S1);

    %% 2. Metropolis法则判断是否接受新解
    [S1,R] = Metropolis(S1,S2,D,T0);  %Metropolis 抽样算法
    
    %% 3. 记录每次迭代过程的最优路线
    if ((count == 1) || (R < Obj(count-1)))
        Obj(count) = R;           %如果当前温度下最优路程小于上一路程，则记录当前路程
    else
        Obj(count) = Obj(count-1);%如果当前温度下最优路程大于上一路程，则记录上一路程
    end
    track(count,:) = S1;
    
    %% 降温
    T0 = q * T0;     
end
toc;

%% VII. 绘制最优路径图 
DrawPath(track(end,:),city,D);

%% VIII. 优化过程迭代图
figure;
plot(1:count,Obj);
grid on;
xlabel('迭代次数');
ylabel('距离');
title('优化过程最短距离变化曲线');
 
%% IX. 输出最优解的路线和总距离
S = track(end,:);
OutputPath(S);
Rlength = PathLength(D,S);
disp(['最短距离：',num2str(Rlength)]);

function Length = PathLength(D,S)
%% 计算各个体的路径长度
% 输入：
% D     两两城市之间的距离
% S     个体的轨迹
% 输出：
% length 每个个体的路径长度
n = size(S,2);
Length = 0;
for i = 1:(n - 1)
    Length = Length + D(S(i),S(i + 1));
end
Length = Length + D(S(n),S(1));      % 加上出发点到终点的距离

function OutputPath(S)
%% 输出路径函数
% 输入：S 路径
S = [S,S(1)];
N = length(S);
p = num2str(S(1));  % p是字符串
for i = 2:N
    p = [p,'―>',num2str(S(i))];
end
disp(p);

function S2 = NewAnswer(S1)
%% 输入
%   S1:当前解
%% 输出
%   S2:新解
 
N = length(S1);
S2 = S1;
a = round(rand(1,2)*(N-1)+1);    %产生两个随机位置，用来交换
W = S2(a(1));
S2(a(1)) = S2(a(2));
S2(a(2)) = W;   %得到一个新路线
end

function [S,R] = Metropolis(S1,S2,D,T)
%% 输入
% S1：  当前解
% S2:   新解
% D:    距离矩阵（两两城市的之间的距离）
% T:    当前温度
%% 输出
% S：   下一个当前解
% R：   下一个当前解的路线距离
 
R1 = PathLength(D,S1);  %计算S1路线长度
R2 = PathLength(D,S2);  %计算S2路线长度 
dC = R2 - R1;           %计算能量差
if dC < 0               %如果能力降低，接受新路线，计算路线长度
    S = S2;
    R = R2;
elseif( exp(-dC/T) >= rand)   %以exp(-dC/T)概率接受新路线
    S = S2;
    R = R2;
else                   %不接受新路线
    S = S1;
    R = R1;
end
end

function DrawPath(S,citys,D)
%% 画路径函数
%输入
% S  待画路径   
% citys  各城市坐标位置
% 无输出(无返回值)
Length = PathLength(D,S);  % 调用计算距离函数
figure;
plot([citys(S,1);citys(S(1),1)],[citys(S,2);citys(S(1),2)],'o-'); 
for i = 1:size(citys,1)
    text(citys(i,1),citys(i,2),['   ' num2str(i)]);  % 标出城市序号
end
text(citys(S(1),1),citys(S(1),2),'       起点');    % 标出起点城市
text(citys(S(end),1),citys(S(end),2),'       终点');   % 标出终点城市
grid on;
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['模拟退火算法优化路径(最短距离：' num2str(Length) ')']);

function D = Distance(citys)
%% 计算两两城市之间的距离
% 输入 citys  各城市的位置坐标
% 输出 D  两两城市之间的距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:n
    for j = i+1:n
        D(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));
        D(j,i) = D(i,j);
    end
end

蚁群算法

%% 1.清空环境变量
clear;
clc;

%% 2.导入数据
load citys_data.mat;          %数据集的变量名为city

%% 3.计算城市间相互距离
n=size(city,1);
D=zeros(n,n);
for i=1:n
    for j=1:n
        if i~=j
            D(i,j)=sqrt(sum((city(i,:)-city(j,:)).^2));
        else
            D(i,j)=1e-4;
        end
    end
end

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

%% 5.迭代寻找最佳路径
h = waitbar(0,' Time Loop ');
tic;
while iter<=itermax
    %% 5.1随机产生各个蚂蚁的起点城市
    start=zeros(m,1);
    for i=1:m
        temp=randperm(n);
        start(i,:)=temp(1);
    end
    Table(:,1)=start;
    citys_index=1:n;

    %% 5.2产生解空间（路径表）
    %逐个蚂蚁路径选择
    for i=1:m
        %逐个城市路径选择
        for j=2:n
            tabu = Table(i,:); %已访问城市
            allow_index=~ismember(citys_index,tabu);
            allow=citys_index(allow_index);
            P=allow;

            %% 5.3计算城市间转移概率
            for k=1:length(allow) 
                P(k)=Tau(Table(i,j-1),allow(k))^alpha * Eta(Table(i,j-1),allow(k))^beta;
            end

            P=P/sum(P);

            %% 5.4轮盘赌法选择下一个访问城市
            Pc=cumsum(P);
            target_index = find(Pc>=rand);
            target=allow(target_index(1));
            Table(i,j)=target;
        end
    end

    %% 5.5计算各个蚂蚁的路径距离
    Length=zeros(m,1);
    for i=1:m
        for j=1:n-1
            Length(i,1)=Length(i,1)+D(Table(i,j),Table(i,j+1));
        end
        Length(i,1)=Length(i,1)+D(Table(i,end),Table(i,1));
    end

    %% 5.6计算最短路径及平均距离
    if iter==1
        [min_Length,min_index]=min(Length);
        Length_best(iter)=min_Length;
        Route_best(iter,:)=Table(min_index,:);
        Length_ave(iter)=mean(Length);
    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-1)>min_Length
            Route_best(iter,:)=Table(min_index,:);
        else
            Route_best(iter,:)=Route_best(iter-1,:);
        end
    end

    %% 5.7更新信息素
    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);
        end
        Delta_Tau(Table(i,end),Table(i,1))=Delta_Tau(Table(i,end),Table(i,1)) + Q/Length(i);
    end
    Tau=(1-rho)*Tau+Delta_Tau;

    % 5.8迭代次数加1，清空路径表
    iter=iter+1;
    Table=zeros(m,n);
    waitbar(iter / itermax, h, sprintf('%s%%', num2str(100 * iter / itermax)));

end
close(h);
toc;

%% 6.结果显示
best_route=Route_best(end,:);
best_length=Length_best(end,:);
OutputPath(best_route);                     %列出最短路径
disp(['最短距离: ' num2str(best_length)]);  %列出最短距离

%% 7.绘图
figure(1);
plot([city(best_route,1);city(best_route(1),1)],[city(best_route,2);city(best_route(1),2)],'o-');
for i=1:size(city,1)
    text(city(i,1),city(i,2),[' ' num2str(i)]);    % 标出城市序号
end
text(city(best_route(1),1),city(best_route(1),2),'    起点');  % 标出起点城市
text(city(best_route(end),1),city(best_route(end),2),'    终点'); % 标出终点城市
grid on;
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['蚁群算法优化路径(最短距离：' num2str(best_length) ')']);

figure(2);
plot(1:itermax,Length_ave,'r:',1:itermax,Length_best,'b-');
grid on;
legend('平均距离','最短距离');
xlabel('迭代次数');
ylabel('距离');
title('各代最短距离与平均距离对比');

function OutputPath(S)
%% 输出路径函数
% 输入：S 路径
S = [S,S(1)];
N = length(S);
p = num2str(S(1));  % p是字符串
for i = 2:N
    p = [p,'―>',num2str(S(i))];
end
disp(p);