pso_m.m

用标准PSO算法解决电力系统最优化问题

%%%%%%%%%%%%%%%  控制参数  %%%%%%%%%%%%%%%%%%%%
DIMENSION_NUM=40;%定义维数
PIRTICLE_NUM=100;%粒子数目
TIME=50;%循环次数
CIR_TIME=200;%叠代次数
wa0=0.9;
wa1=0.4;
%%%%%%%%%%%%%%%  模型参数  %%%%%%%%%%%%%%%%%%%%
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;%罚函数系数
C1=2;
C2=2;
vmax=(pmax-pmin)/2;
allvaluetemp=0;
%%%多样性参数%%%
xmax=500;
L=2*xmax*sqrt(DIMENSION_NUM);
%%%%%%%%%%%%%%%  程序变量  %%%%%%%%%%%%%%%%%%%%
%作图用的临时变量
pso_result=zeros(TIME,1);
allvaluetemp=0;%保存多次结果的值，以便计算平均值
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);
%%%%%%%%%%%%%%%   程 序   %%%%%%%%%%%%%%%%%%%%
%循环次数
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;
    %初始化局部最优值为上面产生的随机变量，全局最优值为随机地最优值
    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;
    
    % 开始循环
    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,:));
            %对粒子速度进行约束处理
            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,:);
            %对粒子位置进行约束处理
            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; 
      
        %更新每个粒子的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;
    %显示每个节点功率及总功率
    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     %显示平均结果
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;