%-----------用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) len(i,1)=len(i,1)+d(j,k);
L(i,j)=1;
elseif (d(j,2) 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=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)=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)=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) 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) 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)=add(i);
for j=1:3
Y{i,j}=X{i,j};
end
end
if add(i) 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) % 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)