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

clear all; clc; close all
%% 基础数据
bcs=[800 350
    400 350];
?
b=[150.7   140  30
246.3   84  30
263.7   172  20
402.8   120  25
543.4   125  20
644.8   118  30
768.0   85  15
789.7   147  15
905.7   67  15
920.2   126  15
188.4   248  5
276.8   252  10
376.8   248  10
471.0   242  10
559.3   235  10
669.5   229  15
765.1   225  10
836.1   225  10
934.7   195  10
179.7   305  15
231.9   323  7.5
311.6   300  15
389.8   303  7.5
478.2   300  10
576.7   292  10
673.8   285  10
768.0   313  10
872.3   310  15
934.7   245  20
978.1   330  15
195.6   384  12.5
297.1   375  15
394.1   381  15
413.0   345  7.5
556.4   360  2.5
586.9   363  15
681.1   348  15
211.6   461  15
217.4   528  20
311.6   465  10
413.0   455  25
501.4   448  25
588.3   456  5
617.3   430  5
691.2   419  15
778.2   395  15
883.9   376  15
253.6   572  5
346.3   575  25
443.4   555  25
528.9   538  20
628.9   512  15
710.0   509  5
734.7   475  5
912.9   448  30
268.1   656  30
504.3   620  15
652.1   580  10
750.6   565  10
843.4   540  10
949.1   523  10
1043.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.0088
0.9377   -1.0860
0.3103    1.7040
0.0086    0.7412
0.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.0860
1.1734    0.0088
1.1734    0.7412
0.3103    1.7040
0.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,‘.-‘)