qpso_m.m

qpso算法解决电力系统最优化问题

%%%%%%%%%%%%%%%  控制参数  %%%%%%%%%%%%%%%%%%%%
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;