数学建模算法之粒子群算法

数学建模之粒子群算法

粒子群优化算法(PSO)是一种进化计算技术(evolutionary computation)，1995 年由Eberhart 博士和kennedy 博士提出，源于对鸟群捕食的行为研究 。该算法最初是受到飞鸟集群活动的规律性启发，进而利用群体智能建立的一个简化模型。粒子群算法在对动物集群活动行为观察基础上，利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得最优解。
粒子公式：

在找到这两个最优值时，粒子根据如下的公式来更新自己的速度和新的位置：
v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[]) (a)

present[] = present[] + v[] (b)

v[] 是粒子的速度, w是惯性权重,present[] 是当前粒子的位置. pbest[] and gbest[] 如前定义. rand () 是介于（0， 1）之间的随机数. c1, c2 是学习因子. 通常 c1 = c2 = 2.
这里我们主要讲解粒子群算法解决TSP问题
旅行商问题，即TSP问题（Traveling Salesman Problem）又译为旅行推销员问题、货郎担问题，是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市，他必须选择所要走的路径，路径的限制是每个城市只能拜访一次，而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值。

例：根据city.xlsx文件，遍历51座城市找到最短路径最优值

解：算法的实现
粒子的表示：
1.TSP问题的一个解为一个序列，可以表示为一个粒子；
2.速度的表示：用一个序列的交换序列表示粒子的速度。
适应度函数的定义：
3.当前序列的路径长度即为适应度值，通过经纬度坐标计算。
惯性因子的定义：
4.自身的交换序列即惯性因子。
该MATLAB代码部分来自MATLAB智能算法分析
实现如下：

function indiFit=fitness(x,City ,Citydistance)% 该函数用于计算个体适应度值
m=size(x,1);
n=size(City,1);
indiFit=zeros(m,1);
for i=1:mfor j=1:n-1indiFit(i)=indiFit(i)+ Citydistance(x(i,j),x(i,j+1));endindiFit(i)=indiFit(i)+ Citydistance(x(i,1),x(i,n));
endfunction dist=dist(x,D)% 该函数用于计算总距离
n=size(x,2);
dist=0;
for i=1:n-1dist=dist+D(x(i),x(i+1));
end
dist=dist+D(x(1),x(n));clc;clear%导入数据 成3列矩阵，第一列为城市名，第二列为 x 坐标，第三列为 y坐标cityCoor=[data(:,2) data(:,3)];%data为导入矩阵cityCoor城市坐标矩阵figure
plot(cityCoor(:,1),cityCoor(:,2),'ms','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
legend('城市位置')
ylim([4 78])%限定y轴上限为78，下限为4title('城市分布图','fontsize',12)
xlabel('km','fontsize',12)
ylabel('km','fontsize',12)
%ylim([min(cityCoor(:,2))-1 max(cityCoor(:,2))+1])grid on %grid on是打开显示坐标轴网格线n=size(cityCoor,1);            %城市数目
cityDist=zeros(n,n);           %城市距离矩阵
for i=1:nfor j=1:nif i~=jcityDist(i,j)=((cityCoor(i,1)-cityCoor(j,1))^2+...(cityCoor(i,2)-cityCoor(j,2))^2)^0.5;%欧氏距离公式endcityDist(j,i)=cityDist(i,j);end
end
nMax=400;                      %进化次数
indiNumber=4000;               %个体数目
individual=zeros(indiNumber,n);
%^初始化粒子位置
for i=1:indiNumberindividual(i,:)=randperm(n);    
endindiFit=fitness(individual,cityCoor,cityDist);% 计算种群适应度
[value,index]=min(indiFit);
tourPbest=individual;                              %当前个体最优
tourGbest=individual(index,:) ;                    %当前全局最优
recordPbest=inf*ones(1,indiNumber);                %个体最优记录
recordGbest=indiFit(index);                        %群体最优记录
xnew1=individual;L_best=zeros(1,nMax);% 循环寻找最优路径
for N=1:nMax%计算适应度值indiFit=fitness(individual,cityCoor,cityDist);%更新当前最优和历史最优for i=1:indiNumberif indiFit(i)<recordPbest(i)recordPbest(i)=indiFit(i);tourPbest(i,:)=individual(i,:);endif indiFit(i)<recordGbestrecordGbest=indiFit(i);tourGbest=individual(i,:);endend[value,index]=min(recordPbest);recordGbest(N)=recordPbest(index);for i=1:indiNumber% 与个体最优进行交叉c1=unidrnd(n-1); %产生交叉位c2=unidrnd(n-1); %产生交叉位while c1==c2c1=round(rand*(n-2))+1;c2=round(rand*(n-2))+1;endchb1=min(c1,c2);chb2=max(c1,c2);cros=tourPbest(i,chb1:chb2);ncros=size(cros,2);      %删除与交叉区域相同元素for j=1:ncrosfor k=1:nif xnew1(i,k)==cros(j)xnew1(i,k)=0;for t=1:n-ktemp=xnew1(i,k+t-1);xnew1(i,k+t-1)=xnew1(i,k+t);xnew1(i,k+t)=temp;endendendend%插入交叉区域xnew1(i,n-ncros+1:n)=cros;%新路径长度变短则接受dist=0;for j=1:n-1dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));enddist=dist+cityDist(xnew1(i,1),xnew1(i,n));if indiFit(i)>distindividual(i,:)=xnew1(i,:);end% 与全体最优进行交叉c1=round(rand*(n-2))+1;  %产生交叉位c2=round(rand*(n-2))+1;  %产生交叉位while c1==c2c1=round(rand*(n-2))+1;c2=round(rand*(n-2))+1;endchb1=min(c1,c2);chb2=max(c1,c2);cros=tourGbest(chb1:chb2); ncros=size(cros,2);      %删除与交叉区域相同元素for j=1:ncrosfor k=1:nif xnew1(i,k)==cros(j)xnew1(i,k)=0;for t=1:n-ktemp=xnew1(i,k+t-1);xnew1(i,k+t-1)=xnew1(i,k+t);xnew1(i,k+t)=temp;endendendend%插入交叉区域xnew1(i,n-ncros+1:n)=cros;%新路径长度变短则接受dist=0;for j=1:n-1dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));enddist=dist+cityDist(xnew1(i,1),xnew1(i,n));if indiFit(i)>distindividual(i,:)=xnew1(i,:);endc1=round(rand*(n-1))+1;   % 变异操作,产生变异位c2=round(rand*(n-1))+1;   %产生变异位while c1==c2c1=round(rand*(n-2))+1;c2=round(rand*(n-2))+1;endtemp=xnew1(i,c1);xnew1(i,c1)=xnew1(i,c2);xnew1(i,c2)=temp;%新路径长度变短则接受dist=0;for j=1:n-1dist=dist+cityDist(xnew1(i,j),xnew1(i,j+1));enddist=dist+cityDist(xnew1(i,1),xnew1(i,n));if indiFit(i)>distindividual(i,:)=xnew1(i,:);endend[value,index]=min(indiFit);L_best(N)=indiFit(index);tourGbest=individual(index,:); endfigure
plot(L_best)
title('算法训练过程')
xlabel('迭代次数')
ylabel('适应度值')
grid onfigure
hold on
plot([cityCoor(tourGbest(1),1),cityCoor(tourGbest(n),1)],[cityCoor(tourGbest(1),2),...cityCoor(tourGbest(n),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
hold on
for i=2:nplot([cityCoor(tourGbest(i-1),1),cityCoor(tourGbest(i),1)],[cityCoor(tourGbest(i-1),2),...cityCoor(tourGbest(i),2)],'ms-','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')hold on
end
legend('规划路径')
scatter(cityCoor(:,1),cityCoor(:,2));
title('规划路径','fontsize',10)
xlabel('km','fontsize',10)
ylabel('km','fontsize',10)grid on
ylim([4 80])
% % 注意：代码仅供参考，不可直接用于自己的数模论文中
% % 国赛对于论文的查重要求非常严格，代码雷同也算作抄袭

该算法求解的答案不一定为最优解，需多次计算。
目前最短路径为
14 24 43 23 7 26 8 48 6 27 51 46 12 47 18 4 17 37 5 38 49 9 50 16 2 11 32 1 22 31 28 3 36 35 20 29 21 34 30 10 39 33 45 15 44 42 19 40 41 13 25
该路径长度为432.5745.

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！