function [pso F] = pso_2D()
%   FUNCTION PSO  --------USE Particle Swarm Optimization Algorithm
%global present;
% close all;
pop_size = 10;          %   pop_size 種群大小
part_size = 2;          %   part_size 粒子大小,                                                                      ** =n-D
gbest = zeros(1,part_size+1);            %   gbest 當前搜索到的最小的值
max_gen = 80;          %   max_gen 最大迭代次數
region=zeros(part_size,2);  %   設定搜索空間范圍
region=[-3,3;-3,3];          %                                                                                      **每一維設定不同范圍


rand('state',sum(100*clock));   %   重置隨機數發生器狀態
arr_present = ini_pos(pop_size,part_size);   %   present 當前位置,隨機初始化,rand()的范圍為0~1

v=ini_v(pop_size,part_size);             %   初始化當前速度


pbest = zeros(pop_size,part_size+1);      %   pbest 粒子以前搜索到的最優值，最后一列包括這些值的適應度
w_max = 0.9;                            %   w_max 權系數最大值
w_min = 0.4;
v_max = 2;             %                                                                                          **最大速度,為粒子的范圍寬度
c1 = 2;                   %   學習因子
c2 = 2;                   %   學習因子
best_record = zeros(1,max_gen);     %   best_record記錄最好的粒子的適應度。
%  ————————————————————————
%   計算原始種群的適應度,及初始化
%  ————————————————————————
arr_present(:,end)=ini_fit(arr_present,pop_size,part_size);

% for k=1:pop_size
%     present(k,end) = fitness(present(k,1:part_size));  %計算原始種群的適應度
% end

pbest = arr_present;                                        %初始化各個粒子最優值
[best_value best_index] = min(arr_present(:,end));         %初始化全局最優，即適應度為全局最小的值，根據需要也可以選取為最大值
gbest = arr_present(best_index,:);
%v = zeros(pop_size,1);          %   v 速度
%  ————————————————————————
%   迭代
%  ————————————————————————
% global m;
% m = moviein(1000);      %生成幀矩陣
x=[-3:0.01:3];
y=[-3:0.01:3];
z=@(x,y) 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...
    - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...
    - 1/3*exp(-(x+1).^2 - y.^2);
for i=1:max_gen
    grid on;
    %     plot3(x,y,z);
    %     subplot(121),ezmesh(z),hold on,grid on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;
    %     subplot(122),ezmesh(z),view([145,90]),hold on,grid on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;
    ezmesh(z),hold on,grid on,plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;
    drawnow
    F(i)=getframe;

    %     ezmesh(z)
    % %     view([-37,90])
    %     hold on;
    %     grid on;
    %     %   plot(-0.0898,0.7126,'ro');
    %     plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*');                                                  %改為三維
    %    axis([-2*pi,2*pi,-pi,pi,-50,10]);
    %     hold off;
    pause(0.01);
    %     m(:,i) = getframe;        %添加圖形
    w = w_max-(w_max-w_min)*i/max_gen;
    %    fprintf('#  %i 代開始！\n',i);
    %   確定是否對打散已經收斂的粒子群——————————————————————————————
    reset = 0;          %   reset = 1時設置為粒子群過分收斂時將其打散，如果＝1則不打散
    if reset==1
        bit = 1;
        for k=1:part_size
            bit = bit&(range(arr_present(:,k))<0.1);
        end
        if bit==1       %   bit=1時對粒子位置及速度進行隨機重置
            arr_present = ini_pos(pop_size,part_size);   %   present 當前位置,隨機初始化
            v = ini_v(pop_size,part_size);           %   速度初始化
            for k=1:pop_size                                    %   重新計算適應度
                arr_present(k,end) = fitness(arr_present(k,1:part_size));
            end
            warning('粒子過分集中！重新初始化……');      %   給出信息
            display(i);
        end
    end

    for j=1:pop_size
        v(j,:) = w.*v(j,:)+c1.*rand.*(pbest(j,1:part_size)-arr_present(j,1:part_size))...
            +c2.*rand.*(gbest(1:part_size)-arr_present(j,1:part_size));                        %  粒子速度更新 (a)

        %   判斷v的大小，限制v的絕對值小于5————————————————————————————
        c = find(abs(v)>6);                                                                                              %**最大速度設置，粒子的范圍寬度
        v(c) = sign(v(c))*6;   %如果速度大于3.14則，速度為3.14

        arr_present(j,1:part_size) = arr_present(j,1:part_size)+v(j,1:part_size);              %  粒子位置更新 (b)
        arr_present(j,end) = fitness(arr_present(j,1:part_size));

        if (arr_present(j,end)>pbest(j,end))&(Region_in(arr_present(j,:),region))     %   根據條件更新pbest,如果是最小的值為小于號,相反則為大于號
            pbest(j,:) = arr_present(j,:);
        end

    end

    [best best_index] = max(arr_present(:,end));                                      %   如果是最小的值為min,相反則為max

    if best>gbest(end)&(Region_in(arr_present(best_index,:),region))                  %   如果當前最好的結果比以前的好，則更新最優值gbest,如果是最小的值為小于號,相反則為大于號
        gbest = arr_present(best_index,:);
    end

    best_record(i) = gbest(end);

end

pso = gbest;

display(gbest);

% figure;
% plot(best_record);
% movie2avi(F,'pso_2D1.avi','compression','MSVC');


%   ***************************************************************************
%      計算適應度
%   ***************************************************************************
function fit = fitness(present)
fit=3*(1-present(1)).^2.*exp(-(present(1).^2) - (present(2)+1).^2) ...                                          %**需要求極值的函數,本例即peaks函數
    - 10*(present(1)/5 - present(1).^3 - present(2).^5).*exp(-present(1).^2-present(2).^2) ...
    - 1/3*exp(-(present(1)+1).^2 - present(2).^2);


function ini_present=ini_pos(pop_size,part_size)
ini_present = 3*rand(pop_size,part_size+1);        %初始化當前粒子位置,使其隨機的分布在工作空間                         %** 6即為自變量范圍


function ini_velocity=ini_v(pop_size,part_size)
ini_velocity =3/2*(rand(pop_size,part_size));        %初始化當前粒子速度,使其隨機的分布在速度范圍內


function flag=Region_in(pos_present,region)
[m n]=size(pos_present);
flag=1;
for j=1:n-1
    flag=flag&(pos_present(1:j)>=region(j,1))&(pos_present(1:j)<=region(j,2));
end


function arr_fitness=ini_fit(pos_present,pop_size,part_size)
for k=1:pop_size
    arr_fitness(k,1) = fitness(pos_present(k,1:part_size));  %計算原始種群的適應度
end