电动汽车未来大规模发展需要众多公共充电站服务,公共充电站应根据电动汽车分布进行合理布局。给出电动汽车分布的预测方法,采用基于排队论的充电机配置方法,提出公共充电站布 局最优规划的数学模型。采用与充电站布局有相似数学特点的 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]]; endTn=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); %计算适应值endPBest=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 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(P)=FX(P); PBest(:,P)=X(:,P); [Fgbest,r]=min(FPBest); Fm(Iter)=Fgbest; if Fgbest [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');endfor k=1:size(bcs,1) str=num2str(k); text(bcs(k,1)+20,bcs(k,2),str,'FontSize',20,'color','red'); endaxis 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'); endaxis equalvv=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');endaxis([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,'.-')