MATLAB----粒子群算法

%Kmeans Cluster Algorithm Based on Particle Optimization Algorithm
clc;
clear all;
format long;
tic
data=[1702.8 1639.79 2068.74 
1877.93 1860.96 1975.3 
867.81 2334.68 2535.1 
1831.49 1713.11 1604.68 
460.69 3274.77 2172.99 
2374.98 3346.98 975.31 
2271.89 3482.97 946.7
1783.64 1597.99 2261.31 
198.83 3250.45 2445.08
1494.63 2072.59 2550.51
1597.03 1921.52 2126.76
1598.93 1921.08 1623.33 
1243.13 1814.07 3441.07
2336.31 2640.26 1599.63 
354 3300.12 2373.61 
2144.47 2501.62 591.51 
426.31 3105.29 2057.8 
1507.13 1556.89 1954.51 
343.07 3271.72 2036.94 
2201.94 3196.22 935.53
2232.43 3077.87 1298.87 
1580.1 1752.07 2463.04 
1962.4 1594.97 1835.95 
1495.18 1957.44 3498.02 
1125.17 1594.39 2937.73 
24.22 3447.31 2145.01 
1269.07 1910.72 2701.97 
1802.07 1725.81 1966.35 
1817.36 1926.4 2328.79 
1860.45 1782.88 1875.13];
%--------参数设定-----------
N=70;%粒子数
c1=1.6;c2=1.6;%设定学习因子值（速度更新参数）
wmax=0.9;wmin=0.4;%设定惯性权重初始及终止值
M=1600;%最大迭代数
K=4;%类别数
[S D]=size(data);%样本数和特征维数
%--------初始化----------------
for i=1:N
clmat(i,:)=randperm(S);%随机取整数
end
clmat(clmat>K)=fix(rand*K+1);%取整函数
fitt=inf*ones(1,N);%初始化个体最优适应度
fg=inf;%初始化群体最优适应度
fljg=clmat(1,:);%当前最优分类
v=rand(N,K*D);%初始速度
x=zeros(N,K*D);%初始化粒子群位置
y=x;%初始化个体最优解
pg=x(1,:);%初始化群体最优解
cen=zeros(K,D);%类别中心定维
fitt2=fitt;%粒子适应度定维
%------循环优化开始------------
for t=1:M
for i=1:N
   ww = zeros(S,K);%%产生零矩阵
   for ii = 1:S
       ww(ii,clmat(i,ii)) = 1;%加权矩阵，元素非0即1
   end
   ccc=[];tmp=0;
   for j = 1:K
        sumcs = sum(ww(:,j)*ones(1,D).*data);
        countcs = sum(ww(:,j));       
       if countcs==0
          cen(j,:) =zeros(1,D);
       else
         cen(j,:) = sumcs/countcs;  %求类别中心
       end
       ccc=[ccc,cen(j,:)];%串联聚类中心
       aa=find(ww(:,j)==1);
       if length(aa)~=0
            for k=1:length(aa)
              tmp=tmp+(sum((data(aa(k),:)-cen(j,:)).^2));%%适应度计算
            end
       end
   end
   x(i,:)=ccc;
   fitt2(i) = tmp; %%适应度值Fitness value  
end
%更新群体和个体最优解
for i=1:N
        if fitt2(i)<fitt(i) 
            fitt(i)=fitt2(i);
            y(i,:)=x(i,:);%个体最优
            if fitt2(i)<fg
            pg=x(i,:);%群体最优
            fg=fitt2(i);%群体最优适应度
            fljg=clmat(i,:);%当前最优聚类
            end
        end
   end
bfit(t)=fg;%最优适应度记录
w = wmax - t*(wmax-wmin)/M;%更新权重,线性递减权重法的粒子群算法
       for i=1:N  
           %更新粒子速度和位置
            v(i,:)=w*v(i,:)+c1*rand(1,K*D).*(y(i,:)-x(i,:))+c2*rand(1,K*D).*(pg-x(i,:));
            x(i,:)=x(i,:)+v(i,:);
            for k=1:K
            cen(k,:)=x((k-1)*D+1:k*D);%拆分粒子位置，获得K个中心
            end
            %重新归类
            for j=1:S
                    tmp1=zeros(1,K);
                    for k=1:K
                    tmp1(k)=sum((data(j,:)-cen(k,:)).^2);%每个样本关于各类的距离
                    end
                    [tmp2 clmat(i,j)]=min(tmp1);%最近距离归类
            end
      end
end
%------循环结束------------
M     %迭代次数
fljg  %最优聚类输出
fg    %最优适应度输出
figure(1)
plot(bfit);%绘制最优适应度轨迹
xlabel('种群迭代次数');
ylabel('适应度');
title('适应度曲线');
cen  %聚类中心
toc```

```bash

M =

        1600


fljg =

  列 1 至 19

     2     2     4     2     4     3     3     2     4     2     2     2     1     3     4     3     4     2     4

  列 20 至 30

     3     3     2     2     1     1     4     1     2     2     2


fg =

     5.074427169449084e+06


cen =

   1.0e+03 *

   1.583318501485389   1.635223368390102   2.714581718419657
   2.181834667703786   1.404155371279036   2.423318117482111
   2.998901968724053   2.358864196437189   2.063591486634854
   3.502722145507502   2.296065717941481   1.666953189720583

历时 34.520749 秒。```