【优化布局】基于粒子群算法的充电站最优布局

213 阅读4分钟

 电动汽车未来大规模发展需要众多公共充电站服务,公共充电站应根据电动汽车分布进行合理布局。给出电动汽车分布的预测方法,采用基于排队论的充电机配置方法,提出公共充电站布 局最优规划的数学模型。采用与充电站布局有相似数学特点的 Voronoi图划分充电站服务区域,服务区内电动汽车考虑快充随机性,采用排队论 M/M/s模型,以电动汽车排队等候时间为标准确定充电站规模,建立公共充电站布局最优规划模型,​用粒子群算法求解。

clear all; clc; close all
%% 基础数据
bcs=[800 350    400 350];
​
b=[150.7   140  30246.3   84  30263.7   172  20402.8   120  25543.4   125  20644.8   118  30768.0   85  15789.7   147  15905.7   67  15920.2   126  15188.4   248  5276.8   252  10376.8   248  10471.0   242  10559.3   235  10669.5   229  15765.1   225  10836.1   225  10934.7   195  10179.7   305  15231.9   323  7.5311.6   300  15389.8   303  7.5478.2   300  10576.7   292  10673.8   285  10768.0   313  10872.3   310  15934.7   245  20978.1   330  15195.6   384  12.5297.1   375  15394.1   381  15413.0   345  7.5556.4   360  2.5586.9   363  15681.1   348  15211.6   461  15217.4   528  20311.6   465  10413.0   455  25501.4   448  25588.3   456  5617.3   430  5691.2   419  15778.2   395  15883.9   376  15253.6   572  5346.3   575  25443.4   555  25528.9   538  20628.9   512  15710.0   509  5734.7   475  5912.9   448  30268.1   656  30504.3   620  15652.1   580  10750.6   565  10843.4   540  10949.1   523  101043.3   502  10];
​
na=1500; %电动汽车数量
​
alp=0.1;  
b(:,4)=round(alp.*b(:,3)./sum(b(:,3)).*na);%每个需求点平均负荷
%b(23,4)=37;  
​
​
ns=6;
​
mui=0.6;   
Nchz=round(mui.*sum(b(:,4))./ns);
​
​
bm=1.0e+003*[0.0086,0.0088;1.1734,0.0088;1.1734,0.7412;0.0086,0.7412;0.0086,0.0088];
​
​
BL=sqrt(8.2*1.0e6./((max(bm(:,1))-min(bm(:,1)))*(max(bm(:,2))-min(bm(:,2)))));%BL为图坐标与实际坐标的比例,为固定参数。
​
​
%划分成两个区域干嘛呢??????
Area2=1.0e+003 *[0.0086    0.00880.9377   -1.08600.3103    1.70400.0086    0.74120.0086    0.0088];
Area2=[Area2,2.*ones(size(Area2,1),1)];%size(Area2,1)=5 ones(5,1)=一个5*1数组
​
Area1=1.0e+003 *[0.9377   -1.08601.1734    0.00881.1734    0.74120.3103    1.70400.9377   -1.0860];
Area1=[Area1,1.*ones(size(Area2,1),1)];
​
vv=[Area1;Area2];  %10*3数组 
for k=1:size(bcs,1)%k=1:2
    Ai=find(vv(:,3)==k);%在vv的第三列查找等于k的元素返回索引值
    xx=vv(Ai,1);%横坐标,充电需求点负荷。。。等于k的点 列向量
    yy=vv(Ai,2);
    kk=convhull(xx,yy);%计算凸包,kk是一个列向量
    in=inpolygon(b(:,1),b(:,2),xx(kk),yy(kk));
    b(in,5)=k; 
end
%%
​
​
​
Ep=[];
​
for i=1:size(bcs,1)
    gb=b(b(:,5)==i,:); 
    Ep=[Ep;[sum(gb(:,4)),round(mui.*sum(gb(:,4))./ns),i]]; 
end
​
​
Tn=8;   %最优充电站数量
PopSize=20; %种群数量
MaxIter=400;   %迭代次数
c1s=2.5;   c2s=0.5;
c1e=0.5;   c2e=2.5;
w_start=0.9;   
w_end=0.4;     %惯性权重w取值范围
Iter=1;        
​
​
xmax=max(Area1(:,1)); xmin=min(Area1(:,1));
ymax=max(Area1(:,2)); ymin=min(Area1(:,2));
x = xmin + (xmax-xmin).*rand(Tn,PopSize);
y = ymin + (ymax-ymin).*rand(Tn,PopSize);
​
X=[x;y]; 
V=rand(Tn*2,PopSize);
Vmax=0.1*max((xmax-xmin),(ymax-ymin));  
inAr1=find(b(:,5)==1); 
bb=[b(inAr1,1:2),b(inAr1,4)]; 
for pk=1:1:PopSize
    [FX(pk),~,~,~,~,~,~,~,~,~]=VorCostCDEV(X(1:Tn,pk),X(Tn+1:end,pk),bb,bcs(1,:),BL); %计算适应值
end
PBest=X;
FPBest=FX;
[Fgbest,r]=min(FX);
Fm(Iter)=Fgbest;
CF=Fgbest;    
Best=X(:,r);  
FBest(Iter)=Fgbest;
FgNum=0;
while (Iter<=MaxIter)%粒子群算法
    Iter=Iter+1 %迭代次数
    w_now=((w_start-w_end)*(MaxIter-Iter)/MaxIter)+w_end;%惯性权重
        A=repmat(X(:,r),1,PopSize);    
    
    R1=rand(Tn*2,PopSize);
    R2=rand(Tn*2,PopSize);   
   
    c1=c1e+(c1s-c1e)*(1-acos(-2*Iter/(MaxIter+1)+1)/pi);
    c2=c2e+(c2s-c2e)*(1-acos(-2*Iter/(MaxIter+1)+1)/pi);
        
    
    V=w_now*V+c1*R1.*(PBest-X)+c2*R2.*(A-X); %粒子速度更新公式
    
   
    changeRows=V>Vmax;
    V(changeRows)=Vmax;
   
    changeRows=V<-Vmax;
    V(changeRows)=-Vmax;
    
    
    X=X+1.0*V;
​
  
   for pk=1:1:PopSize
       [FX(pk),~,~,~,~,~,~,~,~,~]=VorCostCDEV(X(1:Tn,pk),X(Tn+1:end,pk),bb,bcs(1,:),BL);
   end
​
    
    P=FX<FPBest;    
    FPBest(P)=FX(P);  
    PBest(:,P)=X(:,P); 
    
    [Fgbest,r]=min(FPBest);
    Fm(Iter)=Fgbest;        
    
    if Fgbest<CF 
        [FBest,gr]=min(FPBest);  
        Best=PBest(:,gr);    
        CF=Fgbest;  
        FgNum=0;
    else
        FgNum=FgNum+1; 
    end
    
    
    if FgNum>10    
        SubX=r;
        while SubX==r || SubX==0
            SubX=round(rand*(PopSize));
        end
        r=SubX;
        FgNum=0;
    end
    
end 
FBest  
Best  
Fm' 
x =Best(1:Tn);
y =Best(Tn+1:end);
[Fcost,CCS,fcs,ucs,NchI,Epc,CVT,CVTI,CDL,CDLI]=VorCostCDEV(x,y,bb,bcs(1,:),BL)
​
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%作图
figure(1)
a=imread('map1.png');
imshow(a);   
hold on
[vxT,vyT]=VoronoiT(bcs(:,1),bcs(:,2),0);   
plot(bcs(:,1),bcs(:,2),'ks','linewidth',12);
% plot(vxT,vyT,'k-','linewidth',3);   
​
plot(b(:,1),b(:,2),'k*','linewidth',1)   
plot(bm(:,1),bm(:,2),'k-','linewidth',1)  
​
for k=1:length(b(:,4))
    str=num2str(b(k,4));
    text(b(k,1)-15,b(k,2)+25,str,'FontSize',5,'color','black');
end
​
for k=1:size(bcs,1)
    str=num2str(k);
    text(bcs(k,1)+20,bcs(k,2),str,'FontSize',20,'color','red');
 end
axis equal
​
​
[vx,vy]=voronoi(x,y);
plot(x,y,'b^',vx,vy,'b-','linewidth',3); 
for k=1:length(x)
    str=num2str(k);
    text(x(k),y(k),str,'FontSize',20,'color','red');
 end
axis equal
vv=VoronoiArea([x,y],3);
​
bax=b(:,1:2);
for k=1:length(x)
    Ai=find(vv(:,3)==k);
    xx=vv(Ai,1);
    yy=vv(Ai,2);
    kk=convhull(xx,yy); 
    in=inpolygon(bax(:,1),bax(:,2),xx(kk),yy(kk));
    str=num2str(k);
    text(bax(in,1),bax(in,2),str,'FontSize',18,'color','blue');
end
axis([min(bm(:,1))-3 max(bm(:,1))+3 min(bm(:,2))-3 max(bm(:,2))+3])
​
figure(2)
Iter_t=1:1:MaxIter+1;
plot(Iter_t,Fm,'.-')