The future large-scale development of electric vehicles requires the services of numerous public charging stations, which should be rationally arranged according to the distribution of electric vehicles. The prediction method of electric vehicle distribution is given, the charging device configuration method based on queuing theory is adopted, and the mathematical model of optimal planning of public charging station is proposed. Mathematical characteristics similar to charging station layout has a Voronoi diagram of charging station service area, service in the electric car considering quick charge randomness, using queuing theory model of M/M/s, on the basis of the electric car waiting time determine charging station scale, establish public charging stations layout model of optimal planning, using the particle swarm algorithm.

clear all; clc; Close %% 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; % Number of electric vehicles ALP =0.1; b(:,4)=round(alp.*b(:,3)./sum(b(:,3)).*na); % average load per demand point %b(23,4)=37; ns=6; Mui = 0.6; Nchz=round(mui.*sum(b(:,4))./ns); Bm = 1.0 e+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.0 e6. / ((Max (bm (:, 1)) - min (bm (:, 1))) * (Max (bm (:, 2)) - min (bm (:, 2))))); %BL is the ratio of graph coordinates to actual coordinates, and is a fixed parameter. % divided into two regions why ？？？？？？ 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)= a 5*1 array 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 array for k = 1: size (BCS, 1) = % k = 1:2 Ai find (vv (:, 3) = = k); Xx =vv(Ai,1); xx=vv(Ai,1); % x, charging demand point load... The point column vector yy=vv(Ai,2); kk=convhull(xx,yy); % calculation convex hull, kk is a column vector 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; % Optimal number of charging stations PopSize=20; % Population MaxIter=400; % iteration times c1s=2.5; C2s = 0.5; C1e = 0.5; C2e = 2.5; W_start = 0.9; W_end = 0.4; % inertia weight w value range 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)% pSO Iter=Iter+1 % number of iterations w_now=((w_start-w_end)*(MaxIter)/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); % particle velocity update formula 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) % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % drawing 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,'.-')Copy the code