電動汽車未來大規模發展需要衆多公共充電站服務，公共充電站應根據電動汽車分布進行合理布局。給出電動汽車分布的預測方法，采用基于排隊論的充電機配置方法，提出公共充電站布 局最優規劃的數學模型。采用與充電站布局有相似數學特點的 Voronoi圖劃分充電站服務區域，服務區内電動汽車考慮快充随機性，采用排隊論 Ｍ／Ｍ／ｓ模型，以電動汽車排隊等候時間為标準确定充電站規模，建立公共充電站布局最優規劃模型，用粒子群算法求解。
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,'.-')
往期回顧>>>>>>
【模式識别】Matlab指紋識别【優化求解】A*算法解決三維路徑規劃問題 matlab自動識别銀行卡号【優化求解】模拟退火遺傳實作帶時間窗的車輛路徑規劃問題【數學模組化】Matlab實作SEIR模型【優化求解】基于NSGA-2的求解多目标柔性工廠中的房間排程算法【優化求解】蟻群算法求最優值 輸入"源代碼"擷取完整代碼
