%%%%%%%%%%%%%%%  控制參數  %%%%%%%%%%%%%%%%%%%%
DIMENSION_NUM=40;%定義維數
PIRTICLE_NUM=100;%粒子數目
TIME=50;%循環次數
CIR_TIME=200;%疊代次數
wa0=1.0;
wa1=0.5;
%100個粒子 200次 wa0=1.2 wa1=0.4（40）
%20個粒子 1000次 wa0=1.2 wa1=0.55  wa0=1.2 wa1=0.55 （40）
%%%%%%%%%%%%%%%  模型參數  %%%%%%%%%%%%%%%%%%%%
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);
%%%%%%%%%%%%%%%  程序變量  %%%%%%%%%%%%%%%%%%%%
qpso_result=zeros(TIME,1);%保存每次的最好值
allvaluetemp=0;%保存多次結果的值，以便計算平均值
xtemp=zeros(CIR_TIME,1);
ytemp=zeros(CIR_TIME,1);
% qpso_result_pbest_value=zeros(CIR_TIME,TIME,PIRTICLE_NUM);
% qpso_result_pbest=zeros(CIR_TIME,TIME,PIRTICLE_NUM,DIMENSION_NUM);
% qpso_result_x_value=zeros(CIR_TIME,TIME,PIRTICLE_NUM);
% qpso_result_x=zeros(CIR_TIME,TIME,PIRTICLE_NUM,DIMENSION_NUM);
qpso_result_gbest_value=zeros(CIR_TIME,TIME);
qpso_result_gbest=zeros(CIR_TIME,TIME,DIMENSION_NUM);
div=zeros(CIR_TIME,TIME);

%%%%%%%%%%%%%%%   程 序   %%%%%%%%%%%%%%%%%%%%
%循環次數
for time=1:TIME
    %初始化，
    tic;
    %隨機產生每個節點功率
    p=rand(PIRTICLE_NUM,DIMENSION_NUM);
    for i=1:PIRTICLE_NUM
        for j=1 : DIMENSION_NUM
          p(i,j)=rand(1,1)*(pmax(j)-pmin(j))+pmin(j);
        end
    end;
    %初始化局部最優值為上面產生的隨機變量，全局最優值為隨機地最優值
    pbest=p;
    value=zeros(PIRTICLE_NUM,1);
    pbest_value=zeros(PIRTICLE_NUM,1);
    for i=1 : PIRTICLE_NUM
        value(i)=fvalue(pbest(i,:),t);
    end;
    pbest_value=value;
    
    dtemp=1;
    valuetemp=pbest_value(1);
    for i=2 : PIRTICLE_NUM
        if valuetemp > pbest_value(i)
            dtemp=i;
            valuetemp=pbest_value(i);
        end;
    end;%找到最好的那個位置
    gbest=pbest(int32(dtemp),:);%求出全局最好值
    gbest_value=valuetemp;

    % 開始循環
    for cir_time=1 : CIR_TIME
        mbest=sum(pbest)/PIRTICLE_NUM;    
        belt=(wa0-wa1)*(CIR_TIME-cir_time)/CIR_TIME+wa1;
        %計算本次的功率
        for i=1:PIRTICLE_NUM
            fi=rand(1,DIMENSION_NUM);
            p0=fi.*pbest(i,:)+(1-fi).*gbest;
            u=log(1./rand(1,DIMENSION_NUM));
            r=ceil(0.5+rand(1,DIMENSION_NUM));
            p(i,:)=p0+belt*(-1).^r.*abs(mbest-p(i,:)).*u;
            for j=1:DIMENSION_NUM
%                 randtemp=rand(1,1);
%                 p1=randtemp*pbest(i,j)+(1-randtemp)*gbest(j);
%                 u=log(1/rand(1,1));
%                 if rand(1,1)>=0.5
%                     p(i,j)=p1+belt*abs(mbest(j)-p(i,j))*u;
%                 else
%                     p(i,j)=p1-belt*abs(mbest(j)-p(i,j))*u;
%                 end;
                if p(i,j)>pmax(j)
                    p(i,j)=pmax(j);
                end
                if p(i,j)<pmin(j)
                    p(i,j)=pmin(j);
                end
            end;     
        end;
        %更新每個粒子的pbest
        for i=1:PIRTICLE_NUM
            value(i)=fvalue(p(i,:),t);        
            if value(i)<pbest_value(i)
                pbest_value(i)=value(i);
                pbest(i,:)=p(i,:);
            end;
        end;
        %更新gbest
        dtemp=1;
        valuetemp=pbest_value(1);
        for i=2 : PIRTICLE_NUM
            if valuetemp > pbest_value(i)
                dtemp=i;
                valuetemp=pbest_value(i);
            end;
        end;
        %求出全局最好值
        if gbest_value > valuetemp
            gbest=pbest(int32(dtemp),:);
            gbest_value=valuetemp;
        end;
        
        divtemp=diversity(PIRTICLE_NUM,L,DIMENSION_NUM,p,sum(p)/PIRTICLE_NUM);
        div(cir_time,time)=divtemp;
        
%         qpso_result_pbest(cir_time,time,:,:)=pbest;        
%         qpso_result_pbest_value(cir_time,time,:)=pbest_value;
%         qpso_result_x(cir_time,time,:,:)=p;        
%         qpso_result_x_value(cir_time,time,:)=value;
        qpso_result_gbest(cir_time,time,:)=gbest;
        qpso_result_gbest_value(cir_time,time)=gbest_value;
        
        ytemp(cir_time)=ytemp(cir_time)+gbest_value;
        xtemp(cir_time)=cir_time;        
    end;
    %顯示每個節點功率及總功率
    total_power=0.00;
    for i=1 : DIMENSION_NUM
        total_power=total_power+gbest(i);
    end;
    %each_power=gbest
    %total_power
    toc;
    qpso_result(time)=gbest_value;
    allvaluetemp=allvaluetemp+gbest_value;
end;


qpso_result_mean=allvaluetemp/time     %顯示平均結果
qpso_result_process_x=xtemp;
ytemp=ytemp/time;
qpso_result_process_y=ytemp;
savefile='qpso_result.mat';
save(savefile,'qpso_result','qpso_result_mean','qpso_result_process_y','qpso_result_process_x','div','qpso_result_gbest','qpso_result_gbest_value');
%all the date in process
% savefile='qpso_result_large.mat';
% save(savefile,'qpso_result_pbest','qpso_result_pbest_value','qpso_result_x','qpso_result_x_value');

hold on;
plot(xtemp,ytemp,'black');    
        
clear;