pso.m

改进粒子群聚类,通过粒子群寻优算法完成聚类过程

%-----------用PSO算法--------------
function PSO(x)

%------给定初始化条件----------------------------------------------
c1=1.4;             %学习因子1
c2=1.4;             %学习因子2
w=0.8;              %惯性权重
Nmax=100;            %最大迭代次数
D=2;                  %搜索空间维数（未知数个数）
M=20;                  %初始化群体个体数目
[N,N1]=size(x);                %s搜素空间内点数
eps=0.1;           %设置精度(在已知最小值时候用)
a=0;              %设置位置的范围
b=10;
V=[0,10,0,10];       %控制坐标轴的刻度范围

%------初始化种群的个体(可以在这里限定位置和速度的范围)------------
for i=1:N
    for j=1:D
        x(i,j)=a+(b-a)*rand; %随机初始化位置
    end
end

y=x;

%------------初始化各粒子的聚类中心---------------------
C=zeros(20,3);         %记录聚类中心点 
for i=1:20 
    while (C(i,1)==C(i,2))|(C(i,2)==C(i,3))|(C(i,1)==C(i,3))
        C(i,:)=randint(1,3,[1,40]);
    end
end
disp('----------样本点----------')
C
    
%------初始化种群的个体的速度------------
v=cell(20,3);
for i=1:20
    for j=1:3
        v{i,j}=rand(1,2);
    end
end 

X=cell(20,3);  %X{i,j}存储各个样本中各聚类中心的的位置
Y=cell(20,3);  %Y{i,j}存储各个样本中各聚类中心的局部最优位置
Z=cell(1,3);  %Z{i,j}存储所有样本的聚类中心的全局最优位置

%------先计算各个粒子的适应度，并初始化Pi和Pg----------------------
len=zeros(20,3);
L=zeros(20,40);
for i=1:20    %对粒子循环
    for j=1:N   %对给定的样本点循环
        for k=1:3      %对聚类中心
            d(j,k)=(x(j,1)-x(C(i,k),1))^2+(x(j,2)-x(C(i,k),2))^2;
        end
       if (d(j,1)<=d(j,2))&(d(j,1)<=d(j,3))  
           len(i,1)=len(i,1)+d(j,k);
           L(i,j)=1;      
       elseif (d(j,2)<=d(j,3))&(d(j,2)<d(j,1))
           len(i,2)=len(i,2)+d(j,k);         
           L(i,j)=2;
       else
           len(i,3)=len(i,3)+d(j,k);         
           L(i,j)=3;
       end
    end
    l(i)=len(i,1)+len(i,2)+len(i,3);
    fitness1(i)=l(i);             %第i个粒子的局部最优
    pi(i,:)=C(i,:);               %聚类中心
    for j=1:3
        X{i,j}=x(C(i,j),:);
        Y{i,j}=x(C(i,j),:);
    end
end

F=zeros(50,3);
fitness2=fitness1(1);           %全局最优
for i=1:20
    if fitness1(i)<fitness2
        fitness2=fitness1(i);
        pg=C(i,:);        %全局最优的聚类中心位置
        good=i;   %记录第几个粒子达到最好
        for j=1:3
            Z{j}=x(C(i,j),:);
        end        
    end
end
F(1,:)=len(good,:);
fnum=1;

%------显示最好的聚类划分----------------------
figure
for j=1:N
    if L(good,j)==1
        hold on;
        plot(x(j,1),x(j,2),'.b');
    elseif L(good,j)==2
        hold on;
        plot(x(j,1),x(j,2),'.g');
    else
        hold on;
        plot(x(j,1),x(j,2),'.k');
    end
end

for i=1:3
    hold on;
    plot(x(pg(i),1),x(pg(i),2),'or');
end
title('PSO初始聚类划分结果：');
axis(V)
              
   
%------进入主要循环，按照公式依次迭代，直到满足精度要求------------

bianhua=0;    %记录聚类中心是否变化
for t=1:Nmax
    E=zeros(1,20);
    for i=1:20 
        x=y;
        for j=1:3
            x(C(i,j),:)=X{i,j};
        end
        
        if bianhua==0
            %---------求样本点最小距离-------------------
            MinDLength(i)=200;    %样本点间的最小距离  
            for j=1:N-1
                for k=i+1:N
                    Len=(x(j,1)-x(k,1))^2+(x(j,2)-x(k,2))^2;
                    if Len<MinDLength(i)
                        MinDLength(i)=Len;
                    end
                end
            end
            %----------结束-------------------------------
                    
            %----------求聚类中心之间的最小距离-------------
            MinCLength(i)=200;
            for j=1:2   
                for k=j+1:3
                    Len=(x(C(i,j),1)-x(C(i,k),1))^2+(x(C(i,j),2)-x(C(i,k),2))^2; 
                    if Len<MinCLength(i)
                        MinCLength(i)=Len;
                    end
                end
            end
            %------------结束-------------------------------
            
            
            %--------------PSO迭代------------------------
            if MinCLength(i)>MinDLength(i)
                %------------更新粒子速度和位置-----------
                for j=1:3 
                    v{i,j}=w*v{i,j}+c1*rand*(Y{i,j}-X{i,j})+c2*rand*(Z{j}-X{i,j});
                    X{i,j}=X{i,j}+v{i,j};
                end
            else
                bianhua=1;   %标记聚类中心是否变化
                while (C(i,1)==C(i,2))|(C(i,2)==C(i,3))|(C(i,1)==C(i,3))      %重新随机生成类中心位置
                    C(i,:)=randint(3,1,[1,40])
                end
            end
        end
        
        %-----------对该样本重新进行聚类划分，并计算适应度-----------------
        for j=1:3
            x(C(i,j),:)=X{i,j};
        end 
        len(i,:)=[0,0,0];
        for j=1:N   %对给定的样本点循环
            for k=1:3      %对聚类中心
               d(j,k)=(x(j,1)-x(C(i,k),1))^2+(x(j,2)-x(C(i,k),2))^2;
            end
            if (d(j,1)<=d(j,2))&(d(j,1)<=d(j,3))
                E(i)=sqrt(d(j,1))+E(i);
                len(i,1)=d(j,1)+len(i,1);
                L(i,j)=1;      
            elseif (d(j,2)<=d(j,3))&(d(j,2)<d(j,1))
                len(i,2)=d(j,2)+len(i,2);
                E(i)=sqrt(d(j,2))+E(i);
                L(i,j)=2;
            else
                len(i,3)=d(j,3)+len(i,3);
                L(i,j)=3;
                E(i)=sqrt(d(j,3))+E(i);
            end
        end
        
        add(i)=len(i,1)+len(i,2)+len(i,3);
        if add(i)<fitness1(i)
            fitness1(i)=add(i);
            for j=1:3
                Y{i,j}=X{i,j};
            end
        end
        
        if add(i)<fitness2
            fitness2=add(i);
            fgood=E(i);
            fnum=fnum+1;
            F(fnum,:)=fgood;
            for j=1:3
                Z{j}=X{i,j};
            end
            good=i;  %记录达到全局最优的样本为哪个
            goodL=L(i,:);    
            
            %if (F(fnum,1)-F(fnum-1,1)<eps)&(F(fnum,2)-F(fnum-1,2)<eps)&(F(fnum,3)-F(fnum-1,3)<eps)
               % break;
           % end
        end           
       
        
    end

     Emin(t)=fitness2;
     
    %--------显示最好的聚类划分----------------------
     x=y;
     for j=1:3
         x(C(good,j),:)=Z{j};
     end

     if t==10
         figure
         for j=1:N
             if goodL(j)==1
                 hold on;
                 plot(x(j,1),x(j,2),'.b');
             elseif goodL(j)==2
                 hold on;
                 plot(x(j,1),x(j,2),'.g');
             else goodL(j)==3
                 hold on;
                 plot(x(j,1),x(j,2),'.k');
             end
         end
         for j=1:3
             hold on;
             plot(x(C(good,j),1),x(C(good,j),2),'or');
         end
         title('PSO迭代10次后聚类划分结果：');
         axis(V)
     end
     
     if t==20
         figure
         for j=1:N
             if goodL(j)==1
                 hold on;
                 plot(x(j,1),x(j,2),'.b');
             elseif goodL(j)==2
                 hold on;
                 plot(x(j,1),x(j,2),'.g');
             else goodL(j)==3
                 hold on;
                 plot(x(j,1),x(j,2),'.k');
             end
         end
         for j=1:3
             hold on;
             plot(x(C(good,j),1),x(C(good,j),2),'or');
         end
         title('PSO迭代20次后聚类划分结果：');
         axis(V)
     end
     %-----------------------结束---------------------------
     
end
%--------------结束------------------------------------


%----------显示聚类结果--------------------------
disp('-----------PSO迭代次数-------------：')
t
disp('-----------PSO最终聚类中心位置为-------：')
C(good,:)
disp('-----------PSO聚类中心-------------：')
for i=1:3
    Z{i}
end


 x=y;
 for j=1:3
     x(C(good,j),:)=Z{j};
 end
 figure
 for j=1:N
     if goodL(j)==1
         hold on;
         plot(x(j,1),x(j,2),'.b');
     elseif goodL(j)==2
         hold on;
         plot(x(j,1),x(j,2),'.g');
     else goodL(j)==3
         hold on;
         plot(x(j,1),x(j,2),'.k');
     end
 end
 for j=1:3
     hold on;
     plot(x(C(good,j),1),x(C(good,j),2),'or');
 end
 title('PSO聚类划分最好结果：');
 axis(V);
 
 
m=1:t;
figure;
plot(m,Emin,'-');
title('PSO算法每次平均准则函数结果为：');

disp('-----------PSO聚类结果-------------：')
Emin(t)