function [pso F] = pso_2D()
% FUNCTION PSO  --------USE Particle Swarm Optimization Algorithm
% global present;
% close all;
clc;
clear all;
pop_size = 10;                  %   pop_size 種群大小 ///粒子數量
part_size = 2;                 %   part_size 粒子大小 ///粒子的維數 
gbest = zeros(1,part_size+1);   %   gbest 當前搜索到的最小的值 
max_gen = 200;                  %   max_gen 最大迭代次數

%best=zeros(part_size,pop_size*part_size);%xuan


region=zeros(part_size,2);      % 設定搜索空間范圍->解空間
region=10*[-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3;-3,3]; % 每一維設定不同范圍(稱之為解空間，不是可行域空間）

rand('state',sum(100*clock));   % 重置隨機數發生器狀態

%當前種群的信息矩陣，逐代進化的群體 % 當前位置,隨機初始化
% 一個10*3的隨機的矩陣（初始化所有粒子的所有維數的位置值），其中最后一列為  
arr_present = ini_pos(pop_size,part_size); 
          
% 初始化當前速度 
% 一個10*2的隨機的矩陣（初始化所有粒子的所有維數的速度值） 
v=ini_v(pop_size,part_size);  

%不是當前種群，可看作是一個外部的記憶體，存儲每個粒子歷史最優值（2維數值）：根據適應度更新！
%注意：pbest數組10*3  最后一列保存的是適應度
pbest = zeros(pop_size,part_size+1);    % pbest：粒子以前搜索到的最優值，最后一列包括這些值的適應度

% 1*80 保存每代的最優值
best_record = zeros(part_size+1,max_gen);         % best_record數組：記錄每一代的最好的粒子的適應度

w_max = 0.9;          %  w_max權系數最大值
w_min = 0.2;          %  w_min權系數最小值
v_max = 2;            %  最大速度,為粒子的范圍寬度
c1 = 2;               %  學習因子1
c2 = 2;               %  學習因子2

% ————————————————————————
%   計算原始種群的適應度,及初始化
% ————————————————————————
% 注意：傳入的第一個參數是當前的粒子群體 ，ini_fit函數計算每個粒子的適應度
% arr_present(:,end)是最后一列 ，保存每個粒子的適應值，是這樣的！xuan
arr_present(:,end)= ini_fit( arr_present, pop_size, part_size );

% 數組賦值，初始化每個粒子個體的歷史最優值，以后會更新的 
pbest = arr_present;         % 初始化各個粒子最優值

% 找到當前群體中適應度最小的（在最后一列中尋找），best_value
% 改為max，表示關聯度最大
[best_value best_index] = max(arr_present(:,end)); %初始化全局最優，即適應度為全局最小的值，根據需要也可以選取為最大值

% 唯一的全局最優值，是當前代所有粒子中最好的一個
gbest = arr_present(best_index,:);


% 因為是多目標，因此這個-----------------
% 只是示意性的畫出3維的
%x=[-3:0.01:3];
%y=[-3:0.01:3];
%[X,Y]=meshgrid(x,y);
%Z1=(-10)*exp((-0.2)*sqrt(X^2+Y^2));
%Z2=(abs(X))^0.8+abs(Y)^0.8+5*sin(X^3)+5*sin(Y^3);

%z1=@(x,y)(-10)*exp((-0.2)*sqrt(x^2+y^2));
%z2=@(x,y)(abs(x))^0.8+abs(y)^0.8+5*sin(x^3)+5*sin(y^3);
%ezmeshc(z1);grid on;
%ezmeshc(z2);grid on;



%開始進化，直到最大代數截至
for i=1:max_gen
    %grid on;
    %三維圖象 %多維圖象是畫不出來的
    %ezmesh(z),hold on,grid on;
    %畫出粒子群
    %plot3(arr_present(:,1),arr_present(:,2),arr_present(:,3),'*'),hold off;
    %drawnow 
    %flush
    %pause(0.01);
    
    w = w_max-(w_max-w_min)*i/max_gen; % 線形遞減權重

    % 當前進化代數：對于每個粒子進行更新和評價----->>>>>>>
    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的絕對值小于20———————————————————
        for k=1:part_size
            if abs(v(j,k))>20
                rand('state',sum(100*clock));
                v(j,k)=20*rand();
            end
        end

        %前幾列是位置信息
        arr_present(j,1:part_size) = arr_present(j,1:part_size)+v(j,1:part_size);% 粒子位置更新 (b)
        %最后一列是適應度
        arr_present(j,end) = fitness(part_size,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



    %------------混沌---------------------------------
    xlhd = gbest(1:part_size);    
    if(1)
       for p=1:25 %次數
        %1生成
              cxl=rand(1,part_size);
              for j=1:part_size
                  if cxl(j)==0
                      cxl(j)=0.1;
                  end
                  if cxl(j)==0.25
                      cxl(j)=0.26;
                  end
                  if cxl(j)==0.5
                      cxl(j)=0.51;
                  end
                  if cxl(j)==0.75
                      cxl(j)=0.76;
                  end
                  if cxl(j)==1
                      cxl(j)=0.9;
                  end
              end
        %2映射
              al=-30;bl=30;
              rxl=al+(bl-al)*cxl;
        %3搜索
              bate = 0.1;
              xlhd=xlhd+bate*rxl;
              if fitness(part_size,xlhd)>gbest(end)
                  gbest(1:part_size)=xlhd;
                  gbest(end)=fitness(part_size,xlhd);
              end
        %4更新
             for j=1:part_size
                  cxl(j)=4*cxl(j)*(1-cxl(j));
             end
       end
    end
    %-------------混沌--------------------------------
    
    
    
    
    %當前代的最優粒子的適應度（取自）保存
    best_record(:,i) = gbest; % gbest：一個行向量
end

pso = gbest; % 最優個體
display(gbest);

figure;
plot(best_record(end,:));% 最優解與代數的進化關系圖


best=zeros(part_size,max_gen);
for i=1:part_size-1
    best(i,:)=best_record(i,:);
end
pareto1= zeros(1,max_gen);
pareto2= zeros(1,max_gen);
for i=1:max_gen
    pareto1(i)=f1(part_size, best(:,i) );
    pareto2(i)=f2(part_size, best(:,i) );
end

figure;
i=1:max_gen;
%plot(i,pareto1(i),'r*',i,pareto2(i),'g*');
plot(pareto1(i),pareto2(i),'r+');
xlabel('f1');ylabel('f2');
title('Pareto曲線');

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






%子函數
%-------------------------------------------------------------------------
%-------------------------------------------------------------------------
%返回隨機的位置
function ini_present=ini_pos(pop_size,part_size)
ini_present = 10*3*rand(pop_size,part_size+1);       %初始化當前粒子位置,使其隨機的分布在工作空間 
%返回一個隨機的矩陣，10*（2+1），最后一列將用來保存適應度

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


%判斷是否處于范圍內
function flag = Region_in(pos_present,region)
[m n]=size(pos_present); % 1*11  n返回解的維數10
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(part_size,pos_present(k,1:part_size));  %計算原始種群的適應度    
end

%***************************************************************************
%    計算適應度
%***************************************************************************
function fit = fitness(n,xp)
%需要求極值的函數,本例即peaks函數
%y0=[-85.4974,-29.9217]; % 注意：這是基準序列，也就是單個最優的極值
y0=[-9.9907,-7.7507];
%y0=[-39.6162,-18.4561];
%y0=[-86.8312,-29.9217];
y1=[f1(n,xp),f2(n,xp)];  % n為粒子維數
fit=graydegree(2,y0,y1); % 關聯度在某種意義上就是適應度


%目標函數1
function r=f1(n,x)
r=0;
for i=1:n-1
    r=r+(-10)*exp((-0.2)*sqrt(x(i)^2+x(i+1)^2));
end

%目標函數2
function r=f2(n,x)
r=0;
for i=1:n
    r=r+(abs(x(i)))^0.8+5*sin(x(i)^3);
end

%約束函數1
function r=g1(n,x)
r=0;
for i=1:n
    r=0;
end

%約束函數2
function r=g2(n,x)
r=0;
for i=1:n
    r=0;
end

%  灰色關聯度計算函數 ( 越大相似性越好 ) 
%  tn目標函數個數   x0基準序列（一組值）   x1貸檢（一組值）
function gama = graydegree( tn,y0,y1 )
gama=0;
rou =0.5;
kesa= zeros(tn,1);    
m1= abs(y0(1)-y1(1)) ;
m2= abs(y0(1)-y1(1)) ;
for i=1:tn
    if( abs(y0(i)-y1(i))<m1 )
        m1= abs(y0(i)-y1(i));
    end
    if( abs(y0(i)-y1(i))>m2 )
        m2= abs(y0(i)-y1(i));
    end
end
for i=1:tn
    kesa(i) = ( m1+rou*m2)/( abs(y0(i)-y1(i)) +rou*m2 );
    gama = gama + kesa(i);
end
gama = gama/tn;


% 可行解的判決函數  gn為約束條件的個數(暫時未用)   n為解（粒子）的維數
function bool = feasible( x,n )
r=0;
%for i=1:gn
    r=max( 0, g1(n,x), g2(n,x) );%判斷約束條件
%end
if(r>0)
    bool=0; %不可行解
else
    bool=1; %可行解
end