粒子群算法源程序

二维粒子群matlab源程序%function [pso F] = pso_2D()% FUNCTION PSO --------USE Particle Swarm Optimization Algorithm% global present;% close all;clc;clear all;pop_size = 10; % pop_size 种群大小 ///粒子数量part_size = 2; % part_size 粒子大小 ///粒子的维数gbest = zeros(1,part_size+1); % gbest 当前搜索到的最小的值max_gen = 200; % max_gen 最大迭代次数%best=zeros(part_size,pop_size*part_size);%xuanregion=zeros(part_size,2); % 设定搜索空间范围->解空间region=10*[-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3]; % 每一维设定不同范围(称之为解空间，不是可行域空间）rand('state',sum(100*clock)); % 重置随机数发生器状态%当前种群的信息矩阵，逐代进化的群体 % 当前位置,随机初始化% 一个10*3的随机的矩阵（初始化所有粒子的所有维数的位置值），其中最后一列为arr_present = ini_pos(pop_size,part_size);% 初始化当前速度% 一个10*2的随机的矩阵（初始化所有粒子的所有维数的速度值）v=ini_v(pop_size,part_size);%不是当前种群，可看作是一个外部的记忆体，存储每个粒子历史最优值（2维数值）：根据适应度更新！%注意：pbest数组10*3 最后一列保存的是适应度pbest = zeros(pop_size,part_size+1); % pbest：粒子以前搜索到的最优值，最后一列包括这些值的适应度% 1*80 保存每代的最优值best_record = zeros(part_size+1,max_gen); % best_record 数组：记录每一代的最好的粒子的适应度w_max = 0.9; % w_max权系数最大值w_min = 0.2; % w_min权系数最小值v_max = 2; % 最大速度,为粒子的范围宽度c1 = 2; % 学习因子1c2 = 2; % 学习因子2% ————————————————————————% 计算原始种群的适应度,及初始化% ————————————————————————% 注意：传入的第一个参数是当前的粒子群体，ini_fit函数计算每个粒子的适应度% arr_present(:,end)是最后一列，保存每个粒子的适应值，是这样的！xuan arr_present(:,end)= ini_fit( arr_present, pop_size, part_size );% 数组赋值，初始化每个粒子个体的历史最优值，以后会更新的pbest = arr_present; % 初始化各个粒子最优值% 找到当前群体中适应度最小的（在最后一列中寻找），best_value% 改为max，表示关联度最大[best_value best_index] = max(arr_present(:,end)); %初始化全局最优，即适应度为全局最小的值，根据需要也可以选取为最大值% 唯一的全局最优值，是当前代所有粒子中最好的一个gbest = arr_present(best_index,:);% 因为是多目标，因此这个-----------------% 只是示意性的画出3维的%x=[-3:0.01:3];%y=[-3:0.01:3];%[X,Y]=meshgrid(x,y);%Z1=(-10)*exp((-0.2)*sqrt(X^2+Y^2));%Z2=(abs(X))^0.8+abs(Y)^0.8+5*sin(X^3)+5*sin(Y^3);%z1=@(x,y)(-10)*exp((-0.2)*sqrt(x^2+y^2));%z2=@(x,y)(abs(x))^0.8+abs(y)^0.8+5*sin(x^3)+5*sin(y^3); %ezmeshc(z1);grid on;%ezmeshc(z2);grid on;%开始进化，直到最大代数截至for i=1:max_gen%grid on;%三维图象 %多维图象是画不出来的%ezmesh(z),hold on,grid on;%画出粒子群%plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),h old off;%drawnow%flush%pause(0.01);w = w_max-(w_max-w_min)*i/max_gen; % 线形递减权重% 当前进化代数：对于每个粒子进行更新和评价----->>>>>>>for j=1:pop_sizev(j,:) =w.*v(j,:)+c1.*rand.*(pbest(j,1:part_size)-arr_present(j,1:part_size ))...+c2.*rand.*(gbest(1:part_size)-arr_present(j,1:part_size )); % 粒子速度更新 (a)% 判断v的大小，限制v的绝对值小于20———————————————————for k=1:part_sizeif abs(v(j,k))>20rand('state',sum(100*clock));v(j,k)=20*rand();endend%前几列是位置信息arr_present(j,1:part_size) =arr_present(j,1:part_size)+v(j,1:part_size);% 粒子位置更新(b)%最后一列是适应度arr_present(j,end) =fitness(part_size,arr_present(j,1:part_size)); % 适应度更新（保存至最后一列）% 适应度评价与可行域限制if(arr_present(j,end)>pbest(j,end))&(Region_in(arr_present(j,:),regi on)) % 根据条件更新pbest,如果是最小的值为小于号,相反则为大于号 pbest(j,:) = arr_present(j,:); % 更新个体的历史极值 endend% 以下更新全局的极值[best best_index] = max(arr_present(:,end)); % 如果是最小的值为min,相反则为maxif best>gbest(end) &( Region_in(arr_present(best_index,:),region) ) % 如果当前最好的结果比以前的好，则更新最优值gbest,如果是最小的值为小于号,相反则为大于号gbest = arr_present(best_index,:); % 全局的极值end%------------混沌---------------------------------xlhd = gbest(1:part_size);if(1)for p=1:25 %次数%1生成cxl=rand(1,part_size);for j=1:part_sizeif cxl(j)==0cxl(j)=0.1;endif cxl(j)==0.25cxl(j)=0.26;endif cxl(j)==0.5cxl(j)=0.51;endif cxl(j)==0.75cxl(j)=0.76;endif cxl(j)==1cxl(j)=0.9;endend%2映射al=-30;bl=30;rxl=al+(bl-al)*cxl;%3搜索bate = 0.1;xlhd=xlhd+bate*rxl;if fitness(part_size,xlhd)>gbest(end)gbest(1:part_size)=xlhd;gbest(end)=fitness(part_size,xlhd);end%4更新for j=1:part_sizecxl(j)=4*cxl(j)*(1-cxl(j));endendend%-------------混沌--------------------------------%当前代的最优粒子的适应度（取自）保存best_record(:,i) = gbest; % gbest：一个行向量endpso = gbest; % 最优个体display(gbest);figure;plot(best_record(end,:));% 最优解与代数的进化关系图best=zeros(part_size,max_gen);for i=1:part_size-1best(i,:)=best_record(i,:);endpareto1= zeros(1,max_gen);pareto2= zeros(1,max_gen);for i=1:max_genpareto1(i)=f1(part_size, best(:,i) );pareto2(i)=f2(part_size, best(:,i) );endfigure;i=1:max_gen;%plot(i,pareto1(i),'r*',i,pareto2(i),'g*');plot(pareto1(i),pareto2(i),'r+');xlabel('f1');ylabel('f2');title('Pareto曲线');%figure;%plot(,f2(best_record),);% movie2avi(F,'pso_2D1.avi','compression','MSVC');%子函数%------------------------------------------------------------------------- %------------------------------------------------------------------------- %返回随机的位置function ini_present=ini_pos(pop_size,part_size)ini_present = 10*3*rand(pop_size,part_size+1); %初始化当前粒子位置,使其随机的分布在工作空间%返回一个随机的矩阵，10*（2+1），最后一列将用来保存适应度%返回随机的速度function ini_velocity=ini_v(pop_size,part_size)ini_velocity =20*(rand(pop_size,part_size)); %初始化当前粒子速度,使其随机的分布在速度范围内%判断是否处于范围内function flag = Region_in(pos_present,region)[m n]=size(pos_present); % 1*11 n返回解的维数10flag=1;for j=1:n-1flag = flag & ( pos_present(1,j)>=region(j,1) ) &( pos_present(1,j)<=region(j,2) );end%初始化适应度function arr_fitness = ini_fit(pos_present,pop_size,part_size)for k=1:pop_sizearr_fitness(k,1) =fitness(part_size,pos_present(k,1:part_size)); %计算原始种群的适应度%**************************************************** ***********************% 计算适应度%**************************************************** ***********************function fit = fitness(n,xp)%需要求极值的函数,本例即peaks函数%y0=[-85.4974,-29.9217]; % 注意：这是基准序列，也就是单个最优的极值y0=[-9.9907,-7.7507];%y0=[-39.6162,-18.4561];% y0=[-86.8312,-29.9217];y1=[f1(n,xp),f2(n,xp)]; % n为粒子维数fit=graydegree(2,y0,y1); % 关联度在某种意义上就是适应度%目标函数1function r=f1(n,x)r=0;for i=1:n-1r=r+(-10)*exp((-0.2)*sqrt(x(i)^2+x(i+1)^2));end%目标函数2function r=f2(n,x)r=0;for i=1:nr=r+(abs(x(i)))^0.8+5*sin(x(i)^3);%约束函数1function r=g1(n,x)r=0;for i=1:nr=0;end%约束函数2function r=g2(n,x)r=0;for i=1:nr=0;end% 灰色关联度计算函数 ( 越大相似性越好 )% tn目标函数个数x0基准序列（一组值）x1贷检（一组值）function gama = graydegree( tn,y0,y1 )gama=0;rou =0.5;kesa= zeros(tn,1);m1= abs(y0(1)-y1(1)) ;m2= abs(y0(1)-y1(1)) ;for i=1:tnif( abs(y0(i)-y1(i))>m2 ) %------------------应该取大于呢还是小于 m2= abs(y0(i)-y1(i));endendfor i=1:tnkesa(i) = ( m1+rou*m2)/( abs(y0(i)-y1(i)) +rou*m2 );gama = gama + kesa(i);endgama = gama/tn;% 可行解的判决函数 gn为约束条件的个数(暂时未用) n为解（粒子）的维数function bool = feasible( x,n )r=0;%for i=1:gnr=max( 0, g1(n,x), g2(n,x) );%判断约束条件%endif(r>0)bool=0; %不可行解elsebool=1; %可行解endPSO粒子群算法解决旅行商问题的MATLAB源码PSO粒子群算法解决旅行商问题的MATLAB源码%粒子群算法求解旅行商问题% By lReijclose all;clear all;PopSize=500;%种群大小CityNum = 14;%城市数OldBestFitness=0;%旧的最优适应度值Iteration=0;%迭代次数MaxIteration =2000;%最大迭代次数IsStop=0;%程序停止标志Num=0;%取得相同适应度值的迭代次数c1=0.5;%认知系数c2=0.7;%社会学习系数w=0.96-Iteration/MaxIteration;%惯性系数,随迭代次数增加而递减%节点坐标node=[16.47 96.10; 16.47 94.44; 20.09 92.54; 22.39 93.37; 25.23 97.24;...22.00 96.05; 20.47 97.02; 17.20 96.29; 16.30 97.38; 14.05 98.12;...16.53 97.38; 21.52 95.59; 19.41 97.13; 20.09 94.55];%初始化各粒子，即产生路径种群Group=ones(CityNum,PopSize);for i=1:PopSizeGroup(:,i)=randperm(CityNum)';endGroup=Arrange(Group);%初始化粒子速度（即交换序）Velocity =zeros(CityNum,PopSize);for i=1:PopSizeVelocity(:,i)=round(rand(1,CityNum)'*CityNum); %round取整end%计算每个城市之间的距离CityBetweenDistance=zeros(CityNum,CityNum);for i=1:CityNumfor j=1:CityNumCityBetweenDistance(i,j)=sqrt((node(i,1)-node(j,1))^2+(node(i,2)-node(j,2))^2);endend%计算每条路径的距离for i=1:PopSizeEachPathDis(i) = PathDistance(Group(:,i)',CityBetweenDistance);endIndivdualBest=Group;%记录各粒子的个体极值点位置,即个体找到的最短路径IndivdualBestFitness=EachPathDis;%记录最佳适应度值,即个体找到的最短路径的长度[GlobalBestFitness,index]=min(EachPathDis);%找出全局最优值和相应序号%初始随机解figure;subplot(2,2,1);PathPlot(node,CityNum,index,IndivdualBest);title('随机解');%寻优while(IsStop == 0) & (Iteration < MaxIteration)%迭代次数递增Iteration = Iteration +1;%更新全局极值点位置,这里指路径for i=1:PopSizeGlobalBest(:,i) = Group(:,index);end%求pij-xij ,pgj-xij交换序，并以概率c1，c2的保留交换序pij_xij=GenerateChangeNums(Group,IndivdualBest);pij_xij=HoldByOdds(pij_xij,c1);pgj_xij=GenerateChangeNums(Group,GlobalBest);pgj_xij=HoldByOdds(pgj_xij,c2);%以概率w保留上一代交换序Velocity=HoldByOdds(Velocity,w);Group = PathExchange(Group,Velocity); %根据交换序进行路径交换Group = PathExchange(Group,pij_xij);Group = PathExchange(Group,pgj_xij);for i = 1:PopSize % 更新各路径总距离EachPathDis(i) = PathDistance(Group(:,i)',CityBetweenDistance);endIsChange = EachPathDis<IndivdualBestFitness;%更新后的距离优于更新前的，记录序号IndivdualBest(:, find(IsChange)) = Group(:, find(IsChange));%更新个体最佳路径IndivdualBestFitness = IndivdualBestFitness.*( ~IsChange) + EachPathDis.*IsChange;%更新个体最佳路径距离[GlobalBestFitness, index] = min(EachPathDis);%更新全局最佳路径,记录相应的序号if GlobalBestFitness==OldBestFitness %比较更新前和更新后的适应度值;Num=Num+1; %相等时记录加一;elseOldBestFitness=GlobalBestFitness;%不相等时更新适应度值，并记录清零;Num=0;endif Num >= 20 %多次迭代的适应度值相近时程序停止IsStop=1;endBestFitness(Iteration) =GlobalBestFitness;%每一代的最优适应度end%最优解subplot(2,2,2);PathPlot(node,CityNum,index,IndivdualBest);title('优化解');%进化曲线subplot(2,2,3);plot((1:Iteration),BestFitness(1:Iteration));grid on;title('进化曲线');%最小路径值GlobalBestFitnessfunction Group=Arrange(Group)[x y]=size(Group);[NO1,index]=min(Group',[],2); %找到最小值1for i=1:ypop=Group(:,i);temp1=pop([1: index(i)-1]);temp2=pop([index(i): x]);Group(:,i)=[temp2' temp1']';endfunction ChangeNums=GenerateChangeNums(Group,BestVar);[x y]=size(Group);ChangeNums=zeros(x,y);for i=1:ypop=BestVar(:,i);%从BestVar取出一个顺序pop1=Group(:,i);%从粒子群中取出对应的顺序for j=1:x %从BestVar的顺序中取出一个序号NoFromBestVar=pop(j);for k=1:x %从对应的粒子顺序中取出一个序号NoFromGroup=pop1(k);if (NoFromBestVar==NoFromGroup) && (j~=k) %两序号同且不在同一位置ChangeNums(j,i)=k; %交换子pop1(k)=pop1(j);pop1(j)=NoFromGroup;endendendendfunction Hold=HoldByOdds(Hold,Odds)[x,y]=size(Hold);for i=1:xfor j=1:yif rand>OddsHold(i,j)=0;endendendfunction SumDistance=PathDistance(path,CityBetweenDistance)L=length(path); %path为一个循环的节点顺序SumDistance=0;for i=1:L-1SumDistance=SumDistance+CityBetweenDistance(path(i),path(i+1));endSumDistance=SumDistance+CityBetweenDistance(path(1),path(L)); %加上首尾节点的距离function Group=PathExchange(Group,Index)[x y]=size(Group);for i=1:ya=Index(:,i); %取出其中一组交换序pop=Group(:,i); %取出对应的粒子for j=1:x %取出其中一个交换算子作交换if a(j)~=0pop1=pop(j);pop(j)=pop(a(j));pop(a(j))=pop1;endendGroup(:,i)=pop;endfunction PathPlot(node,CityNum,index,EachBest);for i=1:CityNumNowBest(i,:)=node((EachBest(i,index)),:);endNowBest(CityNum+1,:)=NowBest(1,:);plot(node(:,1),node(:,2),'*');line(NowBest(:,1),NowBest(:,2));grid on;。