%%%%%%%%%%%%%%%  控制參數(shù)  %%%%%%%%%%%%%%%%%%%%
DIMENSION_NUM=40;%定義維數(shù)
PIRTICLE_NUM=100;%粒子數(shù)目
TIME=50;%循環(huán)次數(shù)
CIR_TIME=200;%疊代次數(shù)
wa0=0.9;
wa1=0.4;
%%%%%%%%%%%%%%%  模型參數(shù)  %%%%%%%%%%%%%%%%%%%%
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];
%%%%%%%%%%%%%%%  程序參數(shù)  %%%%%%%%%%%%%%%%%%%%
PENISH_COEFFICIENT=100;%罰函數(shù)基數(shù)
t=1;%罰函數(shù)系數(shù)
C1=2;
C2=2;
vmax=(pmax-pmin)/2;
allvaluetemp=0;
%%%多樣性參數(shù)%%%
xmax=500;
L=2*xmax*sqrt(DIMENSION_NUM);
%%%%%%%%%%%%%%%  程序變量  %%%%%%%%%%%%%%%%%%%%
%作圖用的臨時(shí)變量
pso_result=zeros(TIME,1);
allvaluetemp=0;%保存多次結(jié)果的值，以便計(jì)算平均值
xtemp=zeros(CIR_TIME,1);
ytemp=zeros(CIR_TIME,1);
% pso_result_pbest_value=zeros(CIR_TIME,TIME,PIRTICLE_NUM);
% pso_result_pbest=zeros(CIR_TIME,TIME,PIRTICLE_NUM,DIMENSION_NUM);
% pso_result_x_value=zeros(CIR_TIME,TIME,PIRTICLE_NUM);
% pso_result_x=zeros(CIR_TIME,TIME,PIRTICLE_NUM,DIMENSION_NUM);
pso_result_gbest_value=zeros(CIR_TIME,TIME);
pso_result_gbest=zeros(CIR_TIME,TIME,DIMENSION_NUM);
div=zeros(CIR_TIME,TIME);
%%%%%%%%%%%%%%%   程 序   %%%%%%%%%%%%%%%%%%%%
%循環(huán)次數(shù)
for time=1:TIME
    %初始化，
    tic;
    %隨機(jī)產(chǎn)生每個(gè)節(jié)點(diǎn)功率
    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;
    %初始化局部最優(yōu)值為上面產(chǎn)生的隨機(jī)變量，全局最優(yōu)值為隨機(jī)地最優(yōu)值
    for i=1:PIRTICLE_NUM
        for j=1 : DIMENSION_NUM
          v(i,j)=2*rand(1,1)*vmax(j)-vmax(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;
    
    % 開始循環(huán)
    for cir_time=1 : CIR_TIME        
        W=wa0*(CIR_TIME-cir_time)/CIR_TIME+wa1;
        for i=1:PIRTICLE_NUM
%             v(i,:)=W*v(i,:)+C1*rand(1,1)*(pbest(i,:)-p(i,:))+C2*rand(1,1)*(gbest-p(i,:));
            v(i,:)=W*v(i,:)+C1*rand(1,DIMENSION_NUM).*(pbest(i,:)-p(i,:))+C2*rand(1,DIMENSION_NUM).*(gbest-p(i,:));
            %對(duì)粒子速度進(jìn)行約束處理
            for j=1:DIMENSION_NUM
                if v(i,j)>vmax(j)
                    v(i,j)=vmax(j);
                end;
                 if v(i,j)<-vmax(j)
                    v(i,j)=-vmax(j);
                end;          
            end;
            p(i,:)=p(i,:)+v(i,:);
            %對(duì)粒子位置進(jìn)行約束處理
            for j=1:DIMENSION_NUM  
                if p(i,j)>pmax(j)
%                     p(i,j)=p(i,j)-rand(1,1)*(pmax(j)-pmin(j));
                    p(i,j)=pmax(j);
                end
                if p(i,j)<pmin(j)
%                     p(i,j)=p(i,j)+rand(1,1)*(pmax(j)-pmin(j));
                    p(i,j)=pmin(j);
                end
            end             
        end; 
      
        %更新每個(gè)粒子的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;
        
%         pso_result_pbest(cir_time,time,:,:)=pbest;        
%         pso_result_pbest_value(cir_time,time,:)=pbest_value;
%         pso_result_x(cir_time,time,:,:)=p;        
%         pso_result_x_value(cir_time,time,:)=value;
        pso_result_gbest(cir_time,time,:)=gbest;
        pso_result_gbest_value(cir_time,time)=gbest_value;
        
        ytemp(cir_time)=ytemp(cir_time)+gbest_value;
        xtemp(cir_time)=cir_time;       
    end;
    %顯示每個(gè)節(jié)點(diǎn)功率及總功率
    total_power=0.00;
    for i=1 : DIMENSION_NUM
        total_power=total_power+gbest(i);
    end;   
 %   each_power=gbest
 %   total_power
    toc;
    t = toc;
    pso_result(time)=gbest_value;
    allvaluetemp=allvaluetemp+gbest_value;
end;

pso_result_mean=allvaluetemp/time     %顯示平均結(jié)果
pso_result_process_x=xtemp;
ytemp=ytemp/time;
pso_result_process_y=ytemp;
%save data
savefile='pso_result.mat';
save(savefile,'pso_result','pso_result_mean','pso_result_process_y','pso_result_process_x','div','pso_result_gbest','pso_result_gbest_value');
%all the date in process
% savefile='pso_result_large.mat';
% save(savefile,'pso_result_pbest','pso_result_pbest_value','pso_result_x','pso_result_x_value');

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