VRP求解matlab程序

%%theprocedureofantcolonyalgorithmforVRP%%%%%%%%%%%%%initializetheparametersofantcolonyalgorithmscity=[1145215021512641100315926170041302548005128252140061632472100714624640081612428009142239100101632365001114823260012128231120013156217130014129214130015146208300161642089001714120621001814719310001916419390020129189250021155185180022139182700];d=city(:,2:3);x=d(:,1);y=d(:,2);g=city(:,4);m=10;%蚂蚁数mm=max(size(city));alpha=1;belta=3;%决定tao和miu重要性的参数lmda=0;rou=0.4;%衰减系数q0=0.95;%概率tao0=1/(10*400.04);%初始信息素Q=2;%蚂蚁循环一周所释放的信息素defined_phrm=1.0;%initialpheromonelevelvalueQV=6000;%车辆容量vehicle_best=round(sum(g)/QV)+1;%所完成任务所需的最少车数V=40;maxitime=100;best_cost=zeros(1,maxitime);%最好解cost=zeros(maxitime,m);best_tour=[];optcost=inf;opttour=[];maxtt=inf;%计算两点的距离fori=1:mm;forj=1:mm;dist(i,j)=sqrt((d(i,1)-d(j,1))^2+(d(i,2)-d(j,2))^2);end;end;%给taomiu赋初值fori=1:mm;forj=1:mm;ifi~=j;%s(i,j)=dist(i,1)+dist(1,j)-dist(i,j);%节约值tao(i,j)=defined_phrm;%defined_phrm=15.0miu(i,j)=1/dist(i,j);end;end;end;%置deltao(i,j)为0fork=1:mm;fork=1:mm;deltao(i,j)=0;end;end;%%开始迭代forNc=1:maxitime%迭代次数temp=[];fori=1:m%m(10)为蚂蚁数tt=inf;costt=0;sumload=0;cur_pos(i)=1;%从车场出发rn=setdiff(randperm(mm),1);n=1;%A(n)nn=1;part_sol(nn)=1;%部分路径的第一步n_sol=0;%蚂蚁产生的路径数量M_vehicle=500;%是最大车辆数吗t=0;%最佳路径数组的元素数为0hh=1;whilehh~=0%hh为length(rn)whilesumloadQVfork=1:length(rn)ifsumload+g(rn(k))=QV%gama(cur_pos(i),rn(k))=(sumload+g(rn(k)))/QV;%gama越来越大，sumload是个定值A(n)=rn(k);%选择重量适合的点n=n+1;%记录重量适合的点数end;end;na=length(A);ifna==0breakelse%在满足容量的点中计算概率，选择概率大的一点forj=1:nap(j)=10*[(tao(cur_pos(i),A(j)))^alpha]*[(miu(cur_pos(i),A(j)))^belta];%p=compute_prob(A,cur_pos(i),tao,miu,alpha,belta,gama,lmda,i);%？？？endsump=sum(p);p=p/sump;maxp=1e-8;forj=1:naifp(j)maxpmaxp=p(j);index_max=j;end;end;forj=1:naifrand=q0&rand=p(j)index_max=j;breakendendold_pos=cur_pos(i);cur_pos(i)=A(index_max);%当找到一点iftao(old_pos,cur_pos(i))tao(cur_pos(i),old_pos)tao(cur_pos(i),old_pos)=tao(old_pos,cur_pos(i));elsetao(old_pos,cur_pos(i))=tao(cur_pos(i),old_pos);endtao(old_pos,cur_pos(i))=(1-rou)*tao(old_pos,cur_pos(i))+rou*Q;maxtao=10/(1-rou)*Q/200;mintao=maxtao/200;iftao(old_pos,cur_pos(i))maxtaotao(old_pos,cur_pos(i))=maxtao;endiftao(old_pos,cur_pos(i))mintaotao(old_pos,cur_pos(i))=mintao;endnn=nn+1;%nn记录容量和概率都满足的点part_sol(nn)=cur_pos(i);%part_sol为支路sumload=sumload+g(cur_pos(i));temp_load=sumload;%temp_load结束时放出的货物总量rn=setdiff(rn,cur_pos(i));fid=fopen('out_customer.txt','a+');fprintf(fid,'%s%i\t\n','thecurrentpositionis:',cur_pos(i));fprintf(fid,'%s\n','the容量满足的customersetis:');fprintf(fid,'\n\t%i\n',A);fprintf(fid,'------------------------------\n');fclose(fid);na=[];index_max=0;n=1;A=[];endend%whilesumloadQV%当找到一段路%如果当前点为车场n_sol=n_sol+1;%n_sol为路径数,超过5条对其费用加上车辆的派遣费用fid=fopen('out1_solution.txt','a+');fprintf(fid,'%s%i%s','NO.',n_sol,'条路径是:');fprintf(fid,'%i',part_sol);fprintf(fid,'\n');fprintf(fid,'%s','当前送出货物总量是:');fprintf(fid,'%i\n',temp_load);fprintf(fid,'------------------------------\n');fclose(fid);%对所得路径进行路径内3-opt优化（优化函数没给出）fornt=1:length(part_sol);%将所有产生的路径传给一个数组temp(t+nt)=part_sol(nt);%temp为总路径end;t=t+length(part_sol);%总路径元素为tsumload=0;final_sol=setdiff(part_sol,1);rn=setdiff(rn,final_sol);hh=length(rn);part_sol=[];final_sol=[];nn=1;part_sol(nn)=1;A=[];n=1;end%得到一条路的结束点t=t+1;temp(t)=1;ii=length(temp);%计算蚂蚁i产生的路径总长度fori=1:ii-1costt=costt+dist(temp(i),temp(i+1));endcost(Nc,i)=costt;%cost(n_gen,i)=cost_sol(temp,dist)+M_vehicle*(n_sol-vehicle_best);costt=0;ifcost(Nc,i)tttt=cost(Nc,i);best_Nc=Nc;%产生最小费用的代数best_ant=i;%产生最小费用的蚂蚁best_tour=temp;temp=[];end;end%蚂蚁的循环结束点ifmaxttttmaxtt=tt;opttour=best_tour;endbest_cost(1,Nc)=tt;%如果所有蚂蚁均完成一次循环，则用最佳费用所对应的路径对弧进行整体更新%forii=1:mm;%forjj=1:mm;%tao(ii,jj)=(1-rou)*tao(ii,jj)+rou*Q;%end;%end;%对每次迭代的最优蚂蚁进行信息素更新forkk=1:length(best_tour)-1;iftao(best_tour(kk),best_tour(kk+1))tao(best_tour(kk),best_tour(kk+1))tao(best_tour(kk),best_tour(kk+1))=tao(best_tour(kk),best_tour(kk+1));elsetao(best_tour(kk+1),best_tour(kk))=tao(best_tour(kk+1),best_tour(kk));endtao(best_tour(kk),best_tour(kk+1))=tao(best_tour(kk),best_tour(kk+1))+100*Q/tt;end;maxtao=10/(1-rou)*Q/cost(Nc,i);mintao=maxtao/200;fori=1:mmforj=1:mmiftao(i,j)maxtaotao(i,j)=maxtao;endiftao(i,j)mintaotao(i,j)=mintao;endifi==jtao(i,j)=0;endendendfid=fopen('out2_solution.txt','a+');fprintf(fid,'%s%i%s','NO.',n_sol,'路径是:');fprintf(fid,'%i',part_sol);fprintf(fid,'\n');fprintf(fid,'%s%i\n','当前的用户需求量是:',temp_load);fprintf(fid,'%s%f\n','总费用是:',tt);fprintf(fid,'------------------------------\n');fprintf(fid,'\n');fclose(fid);tt=0;end%Nc的结束[optcost,opt_Nc]=min(best_cost)opttourfid=fopen('out3_solution.txt','a+');fprintf(fid,'%s%f\n','总费用是:',optcost);fprintf(fid,'------------------------------\n');fprintf(fid,'%s\n','最终路径是:');fprintf(fid,'%i-',opttour);fprintf(fid,'\n');fclose(fid);%%%作图%%%plot(best_cost);figureplot(x(opttour(1)),y(opttour(1)),'o');holdonfori=1:mmplot(x(i),y(i),'*');holdonendholdonforj=1:length(opttour)-1line([x(opttour(j)),x(opttour(j+1))],[y(opttour(j)),y(opttour(j+1))]);holdonendholdoffxlabel('xcoordinate');ylab