%-------------------------------PSO-----------------------------------%
%                          粒子群優化算法                             %
%---------------------------------------------------------------------%
nn=10;                                    %粒子群規模，即粒子數
m=2;                                      %粒子的維數，即尋優參數的個數
X_pso=rand(nn,m);                         %初始化粒子群的位置
V=rand(nn,m);                             %初始化粒子群的速度
Pbest=rand(nn,m);                              %粒子群的個體最優解
Gbest=rand(1,m);                          %粒子群的全體最優解
%Wmax=0.9;
%Wmin=0.1; 
ITERmax=100;                                      %進化代數
iter=1;
C1=2.05;                                           %學習因子
C2=2.05;                                           %學習因子
%K=2/abs(2-(C1+C2)-sqrt((C1+C2)^2-4*(C1+C2)));      %壓縮因子
%a=1 ;                                              %約束因子
Xmax=100;               
Vmax=10;                                  
while(iter<=ITERmax)
    for i=1:nn
        for j=1:m
            R1=rand;
            R2=rand;
            %W_pso=Wmax-(Wmax-Wmin)/ITERmax*iter;   %慣性因子
            Vinter=V(i,j)+R1*C1*(Pbest(i,j)-X_pso(i,j))+R2*C2*(Gbest(1,j)-X_pso(i,j));
            if Vinter>Vmax
                V(i,j)=Vmax;
            elseif Vinter<-Vmax
                V(i,j)=-Vmax;   
            else 
                V(i,j)=Vinter;
            end
            Xinter=X_pso(i,j)+V(i,j);
            if Xinter>Xmax
                X_pso(i,j)=Xmax;
            elseif Xinter<-Xmax
                X_pso(i,j)=-Xmax; 
            else
                X_pso(i,j)=Xinter;
            end
        end
    end
    for i=1:nn
        pres_X=fit_m(X_pso(i,:));
        pres_Pbest=fit_m(Pbest(i,:));
        pres_Gbest=fit_m(Gbest(1,:));
        if pres_X<pres_Pbest
            Pbest(i,:)=X_pso(i,:);
        end
        if pres_X<pres_Gbest
            Gbest(1,:)=X_pso(i,:)
        end
    end
    err(iter)=pres_Gbest;
    if pres_Gbest<=10^(-30)
        iter=ITERmax;
    else
        iter=iter+1;
    end
end
%figure(1)
ploterr(err/2)
Gbest'
err/2
hold on
% figure(1)
% iter=1:1:ITERmax;
% plot(iter,err(iter),'b-o'); %每次迭代對應的gBest位置

%------------------------------備注-----------------------------------%
%                                                                     %
%---------------------------------------------------------------------%