%%%%%%%%%%%%%%%  控制參數  %%%%%%%%%%%%%%%%%%%%
DIMENSION_NUM=40;%定義維數
PIRTICLE_NUM=100;%粒子數目
TIME=50;%循環次數
CIR_TIME=200;%疊代次數
pc=0.25;
pm=0.1;
%bit_num=20;種群數量,長度
bit_length=21;

%%%%%%%%%%%%%%%  模型參數  %%%%%%%%%%%%%%%%%%%%
P_REQUIRE=8550;%總功
%功率的上下限
pmin=[40  60  80 24 26   68 110 135 135 130    94  94 125 125 125   125 125 220 220 242   242 254 254 254 254   254 254  10  10  10   20 20 20 20 18   18 20 25 25 25 ];
pmax=[80 120 190 42 42  140 300 300 300 300   375 375 500 500 500   500 500 500 500 550   550 550 550 550 550   550 550 150 150 150   70 70 70 70 60   60 60 60 60 60];

%%%%%%%%%%%%%%%  程序參數  %%%%%%%%%%%%%%%%%%%%
PENISH_COEFFICIENT=100;%罰函數基數
t=1;%罰函數系數
%%%多樣性參數%%%
xmax=500;
L=2*xmax*sqrt(DIMENSION_NUM);

%%%%%%%%%%%%%%%  程序變量  %%%%%%%%%%%%%%%%%%%%
ga_result=zeros(TIME,1);%保存每次的最好值
allvaluetemp=0;%保存多次結果的值，以便計算平均值
xtemp=zeros(CIR_TIME,1);
ytemp=zeros(CIR_TIME,1);
% ga_result_x_value=zeros(CIR_TIME,TIME,PIRTICLE_NUM);
% ga_result_x=zeros(CIR_TIME,TIME,DIMENSION_NUM,PIRTICLE_NUM);
ga_result_gbest_value=zeros(CIR_TIME,TIME);
ga_result_gbest=zeros(CIR_TIME,TIME,DIMENSION_NUM);
div=zeros(CIR_TIME,TIME);

for time=1:TIME
    tic;
    %全局最優值
    gbest=zeros(bit_length,DIMENSION_NUM);
    gbestf=zeros(1,DIMENSION_NUM);
    gbest_value=1000000;%假設這個為無窮大正數
    %產生隨機0-1編碼
    x=rand(bit_length,DIMENSION_NUM,PIRTICLE_NUM);
    x=abs(x)>0.5;
    x_float=zeros(DIMENSION_NUM,PIRTICLE_NUM);
    %計算對應的實數值
    f_value=zeros(1,PIRTICLE_NUM);
    for i=1 : PIRTICLE_NUM
        for j=1 : DIMENSION_NUM
            x_float(j,i)=pmin(j)+binary_to_float(x(:,j,i)')*(pmax(j)-pmin(j))/(pow2(bit_length)-1);
        end;
        f_value(i)=fvalue(x_float(:,i)',t);
    end;
    
    %開始循環
    for cirtime=1 : CIR_TIME 
        %計算適應度,結果是累計概率
        fitness=zeros(PIRTICLE_NUM,1);
        fitness=fitness_cal(f_value);
        %進行個體選擇并產生新的種群
        x_new=x;
        for i=1 : PIRTICLE_NUM/2
            selectnum1=rand(1);
            for j=1 : PIRTICLE_NUM
                if selectnum1<fitness(j)
                    selectnum1=j;
                    break;
                end;            
            end;
            selectnum2=rand(1);
            for j=1 : PIRTICLE_NUM
                if selectnum2<fitness(j)
                    selectnum2=j;
                    break;
                end;            
            end;
            %交叉
            x_new=xnum_exchange(x,x_new,selectnum1,selectnum2,i,pc);
            %變異
            x_new=xnum_mutation(x_new,i,pm);
        end;
        x=x_new;
        
        %計算對應的實數值
        for i=1 : PIRTICLE_NUM
            for j=1 : DIMENSION_NUM
                x_float(j,i)=pmin(j)+binary_to_float(x(:,j,i)')*(pmax(j)-pmin(j))/(pow2(bit_length)-1);
            end;
            %計算函數值
            f_value(i)=fvalue(x_float(:,i)',t);
        end;
        
        for i=1 : PIRTICLE_NUM
            if gbest_value>f_value(i)
                %cirtime
                gbest_value=f_value(i);
                gbest=x(:,:,i);
                gbestf=x_float(:,i)';
            end;
        end;
        
        divtemp=diversity(PIRTICLE_NUM,L,DIMENSION_NUM,x_float',sum(x_float')/PIRTICLE_NUM);
        div(cirtime,time)=divtemp;
        
%         ga_result_x(cirtime,time,:,:)=x_float;        
%         ga_result_x_value(cirtime,time,:)=f_value;
        ga_result_gbest(cirtime,time,:)=gbestf;
        ga_result_gbest_value(cirtime,time)=gbest_value;      
        
        ytemp(cirtime)=ytemp(cirtime)+gbest_value;
        xtemp(cirtime)=cirtime;
    end;
    ga_result(time)=gbest_value;
    toc;
end;
allvaluetemp=sum(ga_result);
ga_result_mean=allvaluetemp/TIME;
ga_result_process_x=xtemp;
ytemp=ytemp/time;
ga_result_process_y=ytemp;
%save data
savefile='ga_result.mat';
save(savefile,'ga_result','ga_result_mean','ga_result_process_y','ga_result_process_x','div','ga_result_gbest','ga_result_gbest_value');
%all the date in process
% savefile='ga_result_large.mat';
% save(savefile,'ga_result_x','ga_result_x_value');
% plot the axis
hold on;
plot(xtemp,ytemp,'r');

clear;