de_m.m

微分进化算法来解决电力系统优化问题的源代码

%%%%%%%%%%%%%%%  控制参数  %%%%%%%%%%%%%%%%%%%%
DIMENSION_NUM=40;%定义维数
PIRTICLE_NUM=100;%粒子数目  NP=10*D(DIMENSION_NUM)
TIME=50;%循环次数
CIR_TIME=200;%叠代次数
F=0.5;
CR=0.99;

%%%%%%%%%%%%%%%  模型参数  %%%%%%%%%%%%%%%%%%%%
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;%罚函数系数
N=0;
L=0;
%%%多样性参数%%%
xmax=500;
L1=2*xmax*sqrt(DIMENSION_NUM);

%%%%%%%%%%%%%%%  程序变量  %%%%%%%%%%%%%%%%%%%%
%作图用的临时变量
de_result=zeros(TIME,1);
allvaluetemp=0;%保存多次结果的值，以便计算平均值
xtemp=zeros(CIR_TIME,1);
ytemp=zeros(CIR_TIME,1);
% de_result_x_value=zeros(CIR_TIME,TIME,PIRTICLE_NUM);
% de_result_x=zeros(CIR_TIME,TIME,PIRTICLE_NUM,DIMENSION_NUM);
de_result_gbest_value=zeros(CIR_TIME,TIME);
de_result_gbest=zeros(CIR_TIME,TIME,DIMENSION_NUM);
div=zeros(CIR_TIME,TIME);

%%%%%%%%%%%%%%%   程 序   %%%%%%%%%%%%%%%%%%%%
%循环次数
for time=1:TIME
    %初始化，
    tic;
    %随机产生每个节点功率
    p=zeros(PIRTICLE_NUM,DIMENSION_NUM);
    v=zeros(PIRTICLE_NUM,DIMENSION_NUM);
    vtemp=zeros(PIRTICLE_NUM,DIMENSION_NUM);
    % Initialization
    for i=1:PIRTICLE_NUM
        for j=1 : DIMENSION_NUM
          p(i,j)=rand(1,1)*(pmax(j)-pmin(j))+pmin(j);
        end
    end;
    
    p_value=zeros(PIRTICLE_NUM,1);
    for i=1 : PIRTICLE_NUM
        p_value(i)=fvalue(p(i,:),t);
    end;    
    %求出全局最好的位置
    dtemp=1;
    valuetemp=p_value(1);
    for i=2 : PIRTICLE_NUM
        if valuetemp > p_value(i)
            dtemp=i;
            valuetemp=p_value(i);
        end;
    end;
    gbest=p(int32(dtemp),:);
    gbest_value=valuetemp;
        
    inum=0;
    % 开始循环
    for cir_time=1 : CIR_TIME        
        for inum=1 : PIRTICLE_NUM            
            % Mutation
            randtemp=randint(1,3,PIRTICLE_NUM)+1;
            temp=(sum(randtemp)+inum)/4;
            %检验randremp和inum的数据都不相同
            while ((temp==inum) & (temp==randtemp(1)) & (temp==randtemp(2)) & (temp==randtemp(3))) 
                randtemp=randint(1,3,PIRTICLE_NUM)+1;
                temp=(sum(randtemp)+inum)/4;
            end;
            v(inum,:)=p(randtemp(1),:)+F*(p(randtemp(2),:)-p(randtemp(3),:));
            vtemp(inum,:)=p(inum,:);
            
            % Recombination
            N=randint(1,1,DIMENSION_NUM)+1;
            L=0;
            while (L<=DIMENSION_NUM) & (rand<CR) 
                L=L+1;
            end;
            btemp=N;
            while btemp<(N+L) & (btemp<=DIMENSION_NUM)
                vtemp(inum,btemp)=v(inum,btemp);
                btemp=btemp+1;                
            end;
            
            % Re-calculation
            value=fvalue(vtemp(inum,:),t);
            if value<p_value(inum)
                p_value(inum)=value;
                p(inum,:)=vtemp(inum,:);
                %对粒子位置进行约束处理
                for j=1:DIMENSION_NUM  
                    if p(inum,j)>pmax(j)
%                       p(i,j)=p(i,j)-rand(1,1)*(pmax(j)-pmin(j));
                        p(inum,j)=pmax(j);
                    end
                    if p(inum,j)<pmin(j)
%                       p(i,j)=p(i,j)+rand(1,1)*(pmax(j)-pmin(j));
                        p(inum,j)=pmin(j);
                    end
                end                
            end;
            
        end;

        %求出全局最好的位置
        dtemp=1;
        valuetemp=p_value(1);
        for i=2 : PIRTICLE_NUM
            if valuetemp > p_value(i)
                dtemp=i;
                valuetemp=p_value(i);
            end;
        end;
        gbest=p(int32(dtemp),:);
        gbest_value=valuetemp;
        
        divtemp=diversity(PIRTICLE_NUM,L1,DIMENSION_NUM,p,sum(p)/PIRTICLE_NUM);
        div(cir_time,time)=divtemp;       
        
%         de_result_x(cir_time,time,:,:)=p;        
%         de_result_x_value(cir_time,time,:)=p_value;
        de_result_gbest(cir_time,time,:)=gbest;
        de_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;
    
    de_result(time)=gbest_value;
    allvaluetemp=allvaluetemp+gbest_value;
    
    toc;
end;

de_result_mean=allvaluetemp/time     %显示平均结果
de_result_process_x=xtemp;
ytemp=ytemp/time;
de_result_process_y=ytemp;
savefile='de_result.mat';
save(savefile,'de_result','de_result_mean','de_result_process_y','de_result_process_x','div','de_result_gbest','de_result_gbest_value');
%all the date in process
% savefile='de_result_large.mat';
% save(savefile,'de_result_x','de_result_x_value');
% plot the axis
hold on;
plot(xtemp,ytemp,'y');

clear;