www.gusucode.com > 最优选址的程序 > matlab程序/VorCost_CDEV3.m

    clear all; clc; close all
%% 基础数据
%充电需求点坐标
b =1.0e+003.*[0.1092    0.1348
    0.1197    0.2399
    0.2578    0.1724
    0.4259    0.1739
    0.1257    0.3420
    0.2803    0.3375
    0.4439    0.3360
    0.5505    0.1108
    0.5610    0.2024
    0.5700    0.3075
    0.1332    0.4591
    0.3013    0.4455
    0.4559    0.4380
    0.5850    0.4260
    0.1452    0.6092
    0.3163    0.5341
    0.4739    0.5341
    0.5880    0.5341
    0.3193    0.6407
    0.4784    0.6452
    0.6015    0.6467
    0.6736    0.1574
    0.8657    0.1649
    1.0308    0.1634
    0.6781    0.2924
    0.8327    0.2909
    1.0188    0.2939
    0.6811    0.4425
    0.8192    0.4576
    1.0083    0.4546
    0.6691    0.5431
    0.6976    0.6392
    0.8191    0.6377
    1.0098    0.6347];
%充电需求点常规电力负荷点负荷b(:,3)（kW）
b(:,3)=[2480;2480;8680;11400;890;2340;4160;560;1670;5010;2670;8280;7400;1430;7500;4840;3400;4290;3840;3680;2560;7000;14800;8960;3160;7000;5000;2280;10360;10000;760;6000;7040;5600];

%集中充电站坐标
bcs=[ 937.7296  379.5010
  310.3141  238.4076];



na=4500;

alp=0.1;  
b(:,4)=round(alp.*b(:,3)./sum(b(:,3)).*na);
b(23,4)=37;  


ns=4;

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)))));


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)];

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];   
for k=1:size(bcs,1)
    Ai=find(vv(:,3)==k);
    xx=vv(Ai,1);
    yy=vv(Ai,2);
    kk=convhull(xx,yy);
    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=7;   
PopSize=20; 
MaxIter=300;  
c1s=2.5;   c2s=0.5;
c1e=0.5;   c2e=2.5;
w_start=0.9;   
w_end=0.4;      
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('T3.bmp');
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',5)   
plot(bm(:,1),bm(:,2),'k-','linewidth',3)  

for k=1:length(b(:,4))
    str=num2str(b(k,4));
    text(b(k,1)-15,b(k,2)+25,str,'FontSize',23,'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,'.-')