%function [gBest]=PSO()%PSO主函數(shù)    
pso_size=100;%群體規(guī)模
c1=0.5;c2=0.5;%學習因子
w_max=0.8;%最大權重
w_min=0.4;%最小權重
w=w_max;
Pb=10000000;%適應度值 
Pb1=10000000;
Pb2=10000000;
dimens=1;%待優(yōu)化問題的維數(shù)   
run_max=100;%迭代次數(shù)上限
X=zeros(1,pso_size);
V=zeros(1,pso_size);
Xb=zeros(1,pso_size);
Xb1=zeros(1,pso_size);
Xb2=zeros(1,pso_size);
Yb=zeros(1,pso_size);
Yb1=zeros(1,pso_size);
Yb2=zeros(1,pso_size);
P=zeros(1,8000);
PL=zeros(1,pso_size);
Pp=0;Pp1=0;Pp2=0;
dpBest=zeros(1,pso_size);
dgBest=0;
gBest=0;
s1=zeros(1,run_max);
s2=zeros(1,run_max);
q=0;
L=0.05;%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!精度
t=1;% !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!P(t)中t初值為1
for i=1:pso_size%初始化粒子速度位置
    X(i)=-5+12*rand;
    X1(i)=X(i);
    X2(i)=X(i);
    Xb(i)=X(i);
    Xb1(i)=X(i);
    Xb2(i)=X(i);
    V(i)=(-12+24*rand)/10;
end
for i=1:pso_size   
    Yb1(i)=f1(X(i));%$$$$計算f1中個體適應度值,步驟3$$$$
    if Yb1(i)<Pb1
        Pb1=Yb1(i);%$$$$f1中全局極值$$$$  步驟3,4
        Pp1=X(i);
    end
end
for i=1:pso_size
    Yb2(i)=f2(X(i));%$$$$計算f2中個體適應度值$$$$
    if Yb2(i)<Pb2
        Pb2=Yb2(i);%$$$$f2中全局極值$$$$    步驟3,4
        Pp2=X(i);
    end
end  
gBest=(Pp1+Pp2)/2;%  步驟5
Pp=gBest;
dgBest=abs(Pp1-Pp2);%$$$$計算全局最優(yōu)值的距離
for i=1:pso_size
    dpBest(i)=abs(Xb1(i)-Xb2(i));%$$$$計算各粒子間的距離
    if dpBest(i)<dgBest 
       q=rand;
       if q>0.5
           Xb(i)=Xb1(i);
       else
           Xb(i)=Xb2(i);
       end %$$$$個體極值隨機選取    
    else
        Xb(i)=(Xb1(i)+Xb2(i))/2;%$$$$取兩函數(shù)個體極值的均值
    end
    PL(i)=(Yb1(i)+Yb2(i))/2;%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    if PL(i)<L
    
    P(t)=X(i);
    t=t+1;
    end%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end
%P
for count=1:run_max%循環(huán)迭代
    for k=1:pso_size
            V(k)=w*V(k)+c1*rand*(Xb(k)-X(k))+c2*rand*(Pp-X(k));
            if abs(V(k))>12%判斷是否超出Vmax
                if V(k)>0
                    V(k)=12;
                else
                    V(k)=-12;
                end
            end
             X(k)=X(k)+V(k);
             X1(k)=X(k);
             X2(k)=X(k);
             temp1=f1(X1(k)); %$$$$在f1中計算個體適應度值$$$$
             if temp1<Yb1(k)
                 Yb1(k)=temp1;%$$$$判斷是否更改個體極值$$$$
                 Xb1(k)=X1(k);
             end
             if Yb1(k)<Pb1
                 Pb1=Yb1(k);%$$$$得到新的f1中全局極值Pb1$$$$
                 Pp1=X1(k);
             end
             temp2=f2(X2(k));  %$$$$在f2中計算個體適應度值$$$$
             if temp2<Yb2(k)
                 Yb2(k)=temp2;%$$$$判斷是否更改個體極值$$$$
                 Xb2(k)=X2(k);
             end
             if Yb2(k)<Pb2
                 Pb2=Yb2(k);%$$$$
                 Pp2=X2(k);%$$$$
             end
             gBest=(Pp1+Pp2)/2;%得到全局最優(yōu)
              Pp=gBest;%$$$$全局極值,
             dgBest=abs(Pp1-Pp2);%$$$$計算全局最優(yōu)值的距離
             dpBest(k)=abs(Xb1(k)-Xb2(k));
             if dpBest(k)<dgBest %$$$$比較個體極值的距離$$$$
                 q=rand;
             if q>0.5
              Xb(k)=Xb1(k);
            else
             Xb(k)=Xb2(k);
           end %$$$$個體極值隨機選取  
        else
            Xb(k)=(Xb1(k)+Xb2(k))/2;%$$$$取兩函數(shù)個體極值的均值，得到目標均衡適應度值$$$$
        end
            PL(k)=(Yb1(k)+Yb2(k))/2;
        %flag=0;
        if PL(k)<L%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            a=0;
            for e=1:t-1
                if P(e)==X(k)
                a=1;
                end
            end
            if a~=1
                    P(t)=X(k);
                    t=t+1;
            end
        end
     %   if PL(k)<L        
      %      for e=1:t
       %         if P(t)不等于X(k)
    %P(t)=X(k);
    %t=t+1;
%end
%end
    end%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%end

    w=w_max-count*(w_max-w_min)/run_max;%權值更新公式
    s1(count)=Pb1;
    s2(count)=Pb2;
end
t
for r=1:t
  S3(r)=f1(P(r));
  S4(r)=f2(P(r));
end
s1;
s2;
S3;
S4;
plot(S3,S4,'r*')
hold on
plot(s1,s2,'b#')
hold off