从粒子群算法（PSO）学习鸽群算法(PIO)

鸽群算法和粒子群算法类似都是解决连续空间上的优化问题，而VRP问题都属于组合优化问题，求满足给定约束条件的目标函数最优值的问题，也称离散优化问题。因此需要对鸽群算法公式进行重新定义，得到基于集合的鸽群算法。

鸽群算法：

function Psorout = demo(xy,dmat,Popsize,IterNum)
%利用鸽群优化算法解决TSP问题
nargs = 4;%代表函数要输入参数的个数
for i = nargin:nargs-1
    switch i
        case 0  %产生城市数据
            xy = [488,814;1393,595;2735,2492;4788,4799;4825,1702;789,2927;4853,1120;4786,3757;2427,1276;4002,2530;710,3496;2109,4455;4579,4797;3962,2737;4798,694;3279,747;179,1288;4246,4204;4670,1272;3394,4072;3789,1218;3716,4647;
                1962,1750];
        case 1  %计算距离矩阵
            N = size(xy,1);
            a = meshgrid(1:N);%生成N*N升序矩阵
            dmat = reshape(sqrt(sum((xy(a,:)-xy(a',:)).^2,2)),N,N);% '为矩阵的转置，reshape把数据生成N*N的矩阵
        case 2  %设置粒子群数目
            Popsize = 501;
        case 3  %设置迭代次数
            IterNum = 501;
            IterNum2 = 1200;
        otherwise
    end
end
%对输入参数进行检查
[N,~] = size(xy);%城市个数、维数
[nr,nc] = size(dmat);%距离矩阵的行数和列数
if N~=nr || N~=nc
    error('城市坐标或距离矩阵输入有误')
end

%% PIO参数初始化
R=0.003;                            %磁场参数
pop = zeros(Popsize,N);             %粒子位置
v = zeros(Popsize,N);               %粒子速度
iter = 1;                           %迭代次数计时
iter2 = 1;                          %迭代次数计时
fitness = zeros(Popsize,1);         %适应度函数值
Gbest = zeros(IterNum+IterNum2,N);  %极值路径
Gbest_fitness = zeros(Popsize,1);   %极值
Length_ave = zeros(IterNum,1);
Length_ave2 = zeros(IterNum2,1);
 
%% 产生初始位置和速度
for i = 1:Popsize
    pop(i,:) = randperm(N);
    v(i,:) = randperm(N);
end
%计算粒子适应度值
for i =1:Popsize
    for j =1:N-1
        fitness(i) = fitness(i) + dmat(pop(i,j),pop(i,j+1));
    end
    fitness(i) = fitness(i) + dmat(pop(i,end),pop(i,1));%加上终点回到起点的距离
end
%计算极值
Pbest_fitness = fitness;
Pbest = pop;
[Gbest_fitness(1),min_index] = min(fitness);
Gbest(1,:) = pop(min_index);
Length_ave(1) = mean(fitness);

%% 全局迭代寻优
while iter <IterNum
    %更新迭代次数与惯性系数
    iter = iter +1;
    w = exp(-R*iter); %动态惯性系数
    %更新速度
    %个体极值序列交换部分
    change1 = position_minus_position(repmat(Gbest(iter-1,:),Popsize,1),pop);%将Gbest复制成m行
    change1 = constant_times_velocity(rand,change1);
    %原速度部分
    v = constant_times_velocity(w,v);
    %修正速度
    for i = 1:Popsize
        for j =1:N
            if change1(i,j) ~= 0
                v(i,j) = change1(i,j);
            end
        end
    end
    %更新粒子位置
    pop = position_plus_velocity(pop,v);%更新粒子的位置，也就是更新路径序列
    %适应度值更新
    fitness = zeros(Popsize,1);
    for i = 1:Popsize
        for j = 1:N-1
            fitness(i) = fitness(i) + dmat(pop(i,j),pop(i,j+1));
        end
        fitness(i) = fitness(i) + dmat(pop(i,end),pop(i,1));
    end
    
    for i =1:Popsize
        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(iter-1)
        Gbest_fitness(iter) = minvalue;
        Gbest(iter,:) = pop(min_index,:);
    else
        Gbest_fitness(iter) = Gbest_fitness(iter-1);
        Gbest(iter,:) = Gbest(iter-1,:);
    end
    Length_ave(iter) = mean(fitness);
end
[~,index] = min(Gbest_fitness);
BestRoute = Gbest(index,:);
min(fitness)
figure(1);
plot([xy(BestRoute,1);xy(BestRoute(1),1)],[xy(BestRoute,2);xy(BestRoute(1),2)],'o-')

Np=Popsize;
v2=v;
Pbest2=Pbest;
Pbest_fitness2=Pbest_fitness;
Gbest2=Gbest;
Gbest_fitness2=Gbest_fitness;
%% 局部迭代寻优
while iter2 <IterNum2
    %更新迭代次数
    iter2 = iter2 +1;
    Np=ceil(Np/2);
    sort_fit=sort(fitness);                   %对适应度值进行排序
    fit=sort_fit(Np:length(sort_fit),:);      %选择适应度较好的一半保留
    
    pop2 = zeros(Np,N);
    fitness2 = zeros(Np,1);
    for i=1:Np
        [x,~]=find(fitness==fit(i));
        pop2(i,:)=pop(x(1),:);                                             
        fitness2(i,:)=fitness(x(1),:);  
        Pbest2(i,:)=Pbest(x(1),:);
        Pbest_fitness2(i,:)=Pbest_fitness(x(1),:);
        v2(i,:)=v(x(1),:);
    end
    %更新速度
    change2 = position_minus_position(repmat(Pbest2(1,:),Np,1),pop2);%将Gbest复制成m行
    change2 = constant_times_velocity(rand,change2);

    %修正速度
    for i = 1:Np
        for j =1:N
            if change2(i,j) ~= 0
                v2(i,j) = change2(i,j);
            end
        end
    end
    %更新粒子位置
    pop2 = position_plus_velocity(pop2,v2);%更新粒子的位置，也就是更新路径序列
    %适应度值更新
    fitness2 = zeros(Np,1);
    for i = 1:Np
        for j = 1:N-1
            fitness2(i) = fitness2(i) + dmat(pop2(i,j),pop2(i,j+1));
        end
        fitness2(i) = fitness2(i) + dmat(pop2(i,end),pop2(i,1));
    end
    
    %极值的更新
    for i =1:Np
        if fitness2(i) < Pbest_fitness2(i)
            Pbest_fitness2(i) = fitness2(i);
            Pbest2(i,:) = pop2(i,:);
        end
    end
    [minvalue,min_index] = min(fitness2);
    if minvalue <Gbest_fitness2(IterNum+iter2-2)
        Gbest_fitness2(IterNum+iter2-1) = minvalue;
        Gbest2(IterNum+iter2-1,:) = pop2(min_index,:);
    else
        Gbest_fitness2(IterNum+iter2-1) = Gbest_fitness2(IterNum+iter2-2);
        Gbest2(IterNum+iter2-1,:) = Gbest2(IterNum+iter2-2,:);
    end
    fitness=fitness2;
    pop=pop2;
    Pbest_fitness=Pbest_fitness2;
    Pbest=Pbest2;
    v=v2;
    Length_ave2(iter2-1) = mean(fitness2);
end

[Shortest_Length,index] = min(Gbest_fitness2);
BestRoute = Gbest2(index,:);
min(fitness2)
figure(2);
plot([xy(BestRoute,1);xy(BestRoute(1),1)],[xy(BestRoute,2);xy(BestRoute(1),2)],'o-')
grid on
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(sprintf('鸽群算法优化路径最短距离：%1.2f',Shortest_Length));

目前网上还没有鸽群算法的代码，自己借鉴了网上其他博主关于PSO解决tsp问题的代码，感觉在全局迭代时就已经找到解了，局部迭代好像没啥用，可能是参数设置的问题，反正跑出来的结果不行，请大佬指导！