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