clear all
n = 20;%n維變量
m = 50;%粒子群的數量
c1 = 2;
c2 = 2;
w = 1;
vmax = 0.01;

%隨即產生m個粒子
x = -1+2*rand(m,n);%對每個粒子的每一個變量產生相應的初始位置和速度
v = 0.5*rand(m,n);

%計算適應值

for i = 1:m
    f(i)=obf_pso(x(i,:),n);%每一個粒子的適應度，即為目標函數值。
end
clear i;
%找出個體極值和全局極值
pbx = x;%初值即為個體極值
pbf = f;
[gbf i] = min(pbf);%全局最小值
gbx = pbx(i,:);%
for i =1:m
    v(i,:)=w*v(i,:)+c1*rand*(pbx(i,:)-x(i,:))+c2*rand*(gbx-x(i,:));
    for j=1:n
        if v(i,j)>vmax
            v(i,j)=vmax;
        elseif v(i,j)<vmax
            v(i,j)=-vmax;
        end
    end
    x(i,:) = x(i,:)+v(i,:);
end

%%begin loop
k = 0;
while abs(gbf)>0.001
%while k<20
    for i = 1:m
        f(i)=obf_pso(x(i,:),n);
    end
    %%適應度優于原來的個體則設置當前個體為為最佳個體
    for i=1:m
        if f(i)<pbf(i)
            pbf(i) = f(i);
            pbx(i,:)=x(i,:);
        end
    end
    [gbf i]=min(pbf);
    gbx = pbx(i,:);%找到迭代后全局最佳個體
    for i = 1:m
        v(i,:)=w*v(i,:)+c1*rand*(pbx(i,:)-x(i,:))+c2*rand*(gbx-x(i,:));  %更新速度，并限制在設定的范圍內。
        for j=1:n
            if v(i,j)>vmax
                v(i,j)=vmax;
            elseif v(i,j)<vmax
                v(i,j)=-vmax;
            end
        end
        x(i,:)=x(i,:)+v(i,:);%更新位置
    end
    k = k+1;
    trace(k) = gbf;
    tracea(k,1) = gbx(1);
    tracea(k,2) = gbx(2);
    subplot(1,2,1);
    hold on;
    plot(tracea(k,1),tracea(k,2),'r.','markersize',20);axis([-2,2,-2,2]);
    subplot(1,2,2);
    hold on 
    plot(k,trace(k),'r.')
end
gbx
gbf
k