[优化布局] 基于粒子群算法的充电站优化布局

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,'.-')