www.pudn.com > mtsp-with--pso.zip > ImprovePSO.m, change:2013-07-30,size:8342b


% ===================================================================== 
%                PSO algorithm  
%         Start Date: 2007/10/11 
%       Last Changed: 2007/11/1 
%    Specification: MATLAB 7.01 Programing 
%                           This is PSO algorithm and its application in 
%                           solving the choice of DC 
% 
%线性权重PSO 
%Copyright (c) Zhang Qishan in 2007 
% All rights Reserved 
%========================================================================== 
clear all; 
clc; 
%这里是参数的预备部分,调试时就调整这个里面的参数即可preparation 
nGenNum=2500; %进化代数,就是预设的迭代次数。 
w_fStart=0.95; %0.95;% 粒子开始惯性权重 
w_fEnd=0.4;%粒子结束的惯性权重 
w_fVaryfor=0.7; %粒子惯性权重变化进程次数与总演化次数的比值 
c1=2; %认知部分 加速系数 
c2=2;% 社会部分 加速系数 
c3=1;% 邻域部分 加速系数 
nNhood=3;%邻域宽度 
bNhood=0;%利用邻域模型 
nDisNum=3; % 待选择DC的个数(一定要知道,这是基数,不是序数,即不是需求点的标号) 
nDemandNum=12; % 需求点个数 
rgfDisCenterCapacity=[13 13 13 13 13 13 13 13 13 13 13 13 ];%各配送中心的容量均为13个单位。 
rgfDistoDemDistance=[ 0  1  6  7  4  3  4  6  6  9  8   9; 
    1  0  5  6  5  4  5  7  7 10  9  10; 
    6  5  0  3  6  9 10 12 12 15 14  15; 
    7  6  3  0  3 10 11 13 13 16 15  12; 
    4  5  6  3  0  7  8 10 10 13 12   9; 
    3  4  9 10  7  0  6  4  9 10  6   6; 
    4  5 10 11  8  6  0  2  9  5  4   9; 
    6  7 12 13 10  4  2  0 10  6  2   7; 
    6  7 12 13 10  9  9 10  0  4  8  13; 
    9 10 15 16 13 10  5  6  4  0  4   9; 
    8  9 14 15 12  6  4  2  8  4  0   5; 
    9 10 15 12  9  6  9  7 13  9  5   0];   % rgfDistoDemDistance必须是对称矩阵,表示第i和j点之间的距离。 
%**************************************************************** 
rgfDistoDemUnitCost=[0     1     1     1     1     1     1     1     1     1     1     1 
     1     0     1     1     1     1     1     1     1     1     1     1 
     1     1     0     1     1     1     1     1     1     1     1     1 
     1     1     1     0     1     1     1     1     1     1     1     1 
     1     1     1     1     0     1     1     1     1     1     1     1 
     1     1     1     1     1     0     1     1     1     1     1     1 
     1     1     1     1     1     1     0     1     1     1     1     1 
     1     1     1     1     1     1     1     0     1     1     1     1 
     1     1     1     1     1     1     1     1     0     1     1     1 
     1     1     1     1     1     1     1     1     1     0     1     1 
     1     1     1     1     1     1     1     1     1     1     0     1 
     1     1     1     1     1     1     1     1     1     1     1     0]; % rgfDistoDemUnitCost 必须是对称矩阵,表示第i和j点之间的单位运输成本。 
%**************************************************************** 
rgfDemand=[5 4 2 3 2 4 3 5 4 3 2 2];% rgfdemand需求量 
rgfFixCost=[11 16 14 14 15 13 18 12 11 14 16 11];  %固定费用 
R=1e8;  %罚系数,要求充分大即可 
Particle_size=100; % 例子群的规模 
Vmax=nDisNum-1; 
GM=0; % 最小值目标 
Success=0; %精度达到要求 
ErrGoal =1; 
%********************************************************** 
%下面生成随机初始值x,假设初始的服从[1,3]均匀分布,即[1,m]. 
%初始v 在[-2,2]之间的整数,即[-(m-1),m-1]. 
ub=nDisNum; 
lb=1; 
nInterval=20;% 存储优化的结果 
Iter=0; 
nBestIter=0; 
ParticleS=lb+(ub-lb)*rand(Particle_size,nDemandNum); 
V=round((lb-ub)+2*(ub-lb)*rand(Particle_size,nDemandNum)); 
Particle_demS=lb+(ub-lb)*rand(Particle_size,nDemandNum); 
V_dem=round((lb-ub)+2*(ub-lb)*rand(Particle_size,nDemandNum)); 
%********************************************************** 
%粒子整数规范化 
[Particle_dis,Particle_dem,SelDisNum]=DisDem(ParticleS,Particle_demS,Particle_size,nDemandNum,nDisNum); 
%********************************************************** 
%计算单位运输距离费用矩阵 
rgfDistoDemDisCost=rgfDistoDemDistance.*rgfDistoDemUnitCost;%rgfDistoDemDisCost为单位运输第i至j点的费用矩阵 
%************************************************************************** 
%计算每个粒子的适应度值 
fFEval=zeros(Particle_size,1); 
fFEval=fitness(rgfDistoDemDisCost,rgfFixCost,rgfDemand,Particle_size,SelDisNum,Particle_dem,nDisNum,rgfDisCenterCapacity,R); 
%初始化最佳位置矩阵和相应适应度矩阵 
rgfFPBest=fFEval; 
rgfPDisBestS=ParticleS; 
rgfPDemBestS=Particle_demS; 
%计算粒子群的最优值 
[fFGBest,g]=min(fFEval); 
lastbpf=fFGBest; 
%目前为止粒子群最优配送中心与需求点编号,样本和最优值 
fFBest=fFGBest; 
BestDisS=ParticleS(g,:); 
BestDemS=Particle_demS(g,:); 
BestDis=Particle_dis(g,:); 
BestDem=Particle_dem(g,:); 
history=[0,fFGBest]; 
%生成邻域矩阵 
if bNhood 
    for i=1:Particle_size 
        lo=mod(i-nNhood:i+nNhood,Particle_size); 
        rgnNhood(i,:)=[lo]; 
    end 
    for i=0:nNhood-1 
        rgnNhood(find(rgnNhood==-i))=Particle_size-i; 
    end 
end 
%********************************************************* 
w_fVaryfor=w_fVaryfor*nGenNum; 
InertDec = (w_fStart-w_fEnd)/w_fVaryfor; 
w_fNow=w_fStart; 
tic; 
while ( (Success ==0)&(Iter<=nGenNum)) 
    Iter = Iter+1; 
    % Update the value of the inertia weight w 
    if (Iter<=w_fVaryfor) & (Iter > 1) 
        w_fNow = w_fNow-InertDec; %Change inertia weight 
    end 
    %最佳微粒矩阵 
   A_disS=repmat(ParticleS(g,:),Particle_size,1); 
   A_demS=repmat(Particle_demS(g,:),Particle_size,1); 
   B_DisS=A_disS; 
   B_DemS=A_demS; 
   if bNhood 
       for i=1:Particle_size 
           [fNBest,nb(i)]=min(fFEval(rgnNhood(i,:))); 
           B_DisS(i,:)=ParticleS(rgnNhood(i,nb(i)),:); 
           B_DemS(i,:)=Particle_demS(rgnNhood(i,nb(i)),:); 
       end 
   end 
%生成随机数 
 R1=rand(Particle_size,nDemandNum); 
 R2=rand(Particle_size,nDemandNum); 
  %进化速度 
    if ~bNhood 
        V=w_fNow*V+c1*R1.*(rgfPDisBestS-ParticleS)+c2*R2.*(A_disS-ParticleS); 
        V_dem=w_fNow*V_dem+c1*R1.*(rgfPDemBestS-Particle_demS)+c2*R2.*(A_demS-Particle_demS); 
    else 
        R3=rand(Particle_size,nDemandNum);  
        V=w_fNow*V+c1*R1.*(rgfPDisBestS-ParticleS)+c2*R2.*(A_disS-ParticleS)+c3*R3.*(B_DisS-Particle_demS); 
        V_dem=w_fNow*V_dem+c1*R1.*(rgfPDemBestS-Particle_demS)+c2*R2.*(A_demS-Particle_demS)+c3*R3.*(B_DemS-Particle_demS); 
    end 
    %配送中心粒子更新 
    changeRows=V>Vmax; 
    V(find(changeRows))=Vmax; 
    changeRows=V<-Vmax; 
    V(find(changeRows))=-Vmax; 
    ParticleS=ParticleS+V; 
    changeRows=ParticleS>nDisNum; 
    ParticleS(find(changeRows))=nDisNum; 
    changeRows=ParticleS<1; 
    ParticleS(find(changeRows))=1; 
    %需求点粒子更新 
    changeRows=V_dem>Vmax; 
    V_dem(find(changeRows))=Vmax;; 
    changeRows=V_dem<-Vmax; 
    V_dem(find(changeRows))=-Vmax; 
    Particle_demS=Particle_dem+V; 
    changeRows=Particle_demS>nDisNum; 
    Particle_demS(find(changeRows))=nDisNum; 
    changeRows=Particle_demS<1; 
    Particle_demS(find(changeRows))=1; 
    %配送中心和需求点整数化 
    [Particle_dis,Particle_dem,SelDisNum]=DisDem(ParticleS,Particle_demS,Particle_size,nDemandNum,nDisNum); 
    fFEval=fitness(rgfDistoDemDisCost,rgfFixCost,rgfDemand,Particle_size,SelDisNum,Particle_dem,nDisNum,rgfDisCenterCapacity,R); 
   %更新每个微粒的最佳位置和适应度值 
    changeRows=fFEval<rgfFPBest; 
    rgfFPBest(find(changeRows))=fFEval(find(changeRows)); 
    rgfPDisBestS(find(changeRows),:)=ParticleS(find(changeRows),:); 
    rgfPDemBesS(find(changeRows),:)=Particle_demS(find(changeRows),:); 
    lastbpart_disS=rgfPDisBestS(g,:); 
    lastbpart_demS=rgfPDemBestS(g,:); 
    %求粒子群的最优值并更新索引标识 
    [fFGBest,g]=min(rgfFPBest); 
   %如果适应度改善,则更新最好的 
   if fFGBest<lastbpf 
       [fFBest,b]=min(rgfFPBest); 
       BestDisS=ParticleS(b,:); 
       BestDemS=Particle_demS(b,:); 
       BestDis=Particle_dis(b,:); 
       BestDem=Particle_dem(b,:); 
       nBestIter=Iter; 
   end 
   if nInterval &(rem(Iter,nInterval)==0) 
      history(size(history,1)+1,:)=[Iter,fFBest]; 
   end 
    
    if abs(fFGBest-GM) <= ErrGoal     %GBest 
        Success = 1; 
    elseif abs(fFBest-GM)<=ErrGoal    %Best 
        Success = 1; 
    else 
        lastbpf=fFGBest; 
    end % if  
end %while 
fBestIterTime=toc; 
disp(sprintf('最优迭代时间=%f\t\t  迭代次数=%4d\t\t 最优迭代次数=%4d\t\t 最优值=%f',fBestIterTime, Iter-1,nBestIter,fFGBest)); 
BestDis 
BestDem 
nRows=size(history,1); 
for i=1:nRows 
disp(sprintf('迭代次数=%4d\t\t   最优值=%f\t\t ',history(i,1),history(i,2))); 
end