基于matlab的麻雀算法求解

一麻雀算法讲解

建立麻雀搜索算法的数学模型，主要规则如下所述：

      1.1发现者通常拥有较高的能源储备并且在整个种群中负责搜索到具有丰富食物的区域，为所有的加入者提供觅食的区域和方向。在模型建立中能量储备的高低取决于麻雀个体所对应的适应度值(Fitness Value)的好坏。
     1.2一旦麻雀发现了捕食者，个体开始发出鸣叫作为报警信号。当报警值大于安全值时，发现者会将加入者带到其它安全区域进行觅食。发现者和加入者的身份是动态变化的。只要能够寻找到更好的食物来源，每只麻雀都可以成为发现者，但是发现者和加入者所占整个种群数量的比重是不变的。也就是说，有一只麻雀变成发现者必然有另一只麻雀变成加入者。
     1.3加入者的能量越低，它们在整个种群中所处的觅食位置就越差。一些饥肠辘辘的加入者更有可能飞往其它地方觅食，以获得更多的能量。
     1.4在觅食过程中，加入者总是能够搜索到提供最好食物的发现者，然后从最好的食物中获取食物或者在该发现者周围觅食。与此同时，一些加入者为了增加自己的捕食率可能会不断地监控发现者进而去争夺食物资源。
     1.5当意识到危险时，群体边缘的麻雀会迅速向安全区域移动，以获得更好的位置，位于种群中间的麻雀则会随机走动，以靠近其它麻雀。

二求解结果展示

                                                   14个城市路径求解

三matalb代码展示

%%  基于麻雀算法的tsp问题
%% 1.清空环境变量
clear;
clc;
%% 2.导入数据
% load citys51.mat;          %数据集的变量名为citys
% % %% 3.计算城市间相互距离
% citys=citys51;
load citys_data.mat;   
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
%%4 初始化参数
ST=0.8;    %预警阈值
m=500;       %麻雀种群数量
popn=m;
pop=zeros(m,n);
PD_percent=0.2;   
SD_percent=0.2;      
genmax=1000;       %最大迭代次数
gen=1;
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);     %各代路径的平均长度
%% 5产生随机麻雀种群
% 随机产生麻雀位置
for i=1:m
    pop(i,:)=randperm(n);
end
% 计算适应度函数值
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
% 极端值
% Pbest_fitness=fitness;   %进行预分配
% Pbest=pop;
% [Gbest_fitness(1),min_index]=min(fitness);
% Gbest(1,:)=pop(min_index,:);
PDn=round(PD_percent*popn);    %捕食者个数
SDn=round(SD_percent*popn);    %预警者的个数
while   gen<=genmax
   [~,sortindex]=sort(fitness);  %对麻雀的适应度升序排列，为了更好的找到发现者与跟随者
   fitness_new=fitness; 
   pop_new=pop;
   [fmax,worse]=max(fitness);   %找到最差位置
   pop_worst=pop(worse,:);
   [fmin,best]=min(fitness);    %找到最佳的位置
   pop_best=pop(best,:);
   if rand<ST          %根据预警值判断位置变换
        for i=1:PDn
           a=sortindex(i);
           r=exp(-i/(rand*genmax));
           if rand<=r    %距离近，保持不动
             pop_new(a,:)=pop_new(a,:);
           else
             dim=round(0.2*m);
             pop_new(a,:)=position_change_little(pop(a,:),dim);
           end
        end
    else
        for i=1:PDn
            a=sortindex(i);
            dim=round(rand*m);
            pop_new(a,:)=position_change_little(pop(a,:),dim); 
        end
   end
    fitness_new=fobi_compute(pop_new,D);   %得到新的函数值
    [~,best2]=min(fitness_new);  %找到更新后的值
    pop_best2=pop(best2,:);    %找到最好的位置二
   % 跟随者位置更新
    for i=PDn+1:popn
           a=sortindex(i);
       if i>0.5*popn
           pop_new(a,:)=randperm(n);   %十分饥饿，自寻食物
       else
             change1=position_minus_position(pop_best2,pop_new(a,:));  %记录最好位置到当前位置变换 
             change1=constant_times_velocity(rand,change1);
             pop_new(a,:)=position_plus_velocity(pop_new(a,:),change1);        
       end
    end
    
    fitness_new=fobi_compute(pop_new,D);    %得到新的函数值
    eta=10e-8;
    index=randperm(popn,SDn);  %随机分配预警者 
     
    for i=1:SDn  
        a=index(i);
        if   fitness(a)>fmin
             change1=position_minus_position(pop_best,pop_new(a,:));  %记录最好位置到当前位置变换 
             change1=constant_times_velocity(rand,change1);
             pop_new(a,:)=position_plus_velocity(pop_new(a,:),change1);        %向着中心位置靠拢
        else
             dim=1;
             pop_new(a,:)=position_change_little(pop(a,:),dim);  %中心位置内，小变动
        end   
    end
     fitness_new=fobi_compute(pop_new,D);      %当前函数最佳值
     index=find(fitness_new<fitness);  %两者相互比较，找到最小的
     fitness(index,:)=fitness_new(index);
     pop(index,:)=pop_new(index,:);
     [fmin,best]=min(fitness);    %找到最佳的位置
     Convergence_curve(gen)=fmin;
     result=pop(best,:);
     gen=gen+1;
end

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

function change=position_minus_position(best,pop)
  %记录将pop变成best的交换序列
for i=1:size(best,1)
    for j=1:size(best,2)
        change(i,j)=find(pop(i,:)==best(i,j));
        temp=pop(i,j);
        pop(i,j)=pop(i,change(i,j));
        pop(i,change(i,j))=temp;
    end
end
end

function  change_little=position_change_little(pop,dim)
  %对位置进行晓得变动
  n=size(pop,2);
  for i=1:dim
      a=randperm(n,1);
      b=randperm(n,1);
      tempt=pop(:,a);
      pop(:,a)=pop(:,b);
      pop(:,b)=tempt;  
  end
     change_little=pop;
end

function  fobj=fobi_compute(pop,D)
    m=size(pop,1);  n=size(D);
    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
    fobj=fitness;
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

function pop = position_plus_velocity(pop,v)
%利用速度记录的交换序列进行位置修正
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