A brief introduction to the principle
The basic principles of ant colony algorithm
1. Ants release pheromones along the path.
Encounter hasn’t walked the intersection, randomly choose a road to go. At the same time, pheromones related to path length are released.
3. Pheromone concentration is inversely proportional to path length. When later ants come across the intersection again, they choose the path with higher pheromone concentration.
4. The pheromone concentration on the optimal path increases.
5. Finally, the ant colony finds the optimal feeding path.
There are more ants on the shorter side of the colony than on the longer side, leaving more pheromones behind, and over time the ants only go on the shorter side.
Ant colony algorithm has two main steps to solve TSP :(TSP problem is to find the shortest Hamilton loop)
1. Path construction
Random proportion rule in AS; For each ant K, the path memory vector R^ K records the serial number of all cities k has passed through according to the access order. Suppose ant K’s current city is I, then its probability of selecting city J as the next access object is:
2. Pheromone renewal
Here m is the number of ants, ρ is the evaporation rate of the pheromone, specified 0≤ ρ≤1, in AS usually set to ρ =0.5, δ τij is the amount of pheromone released by the KTH ant on the side it passes, which is equal to the inverse of the length of the construction path of ant K round. C to the k is the length of the path, which is the sum of the lengths of all the edges in R to the k.
Construction diagram: The construction diagram is consistent with the problem description diagram. The set C of components corresponds to the set of points (C=N), and the set of connecting corresponding edges (L=A), and each edge has A weight, representing the distance between points I and j.
Constraints: All cities must be visited and each city must be visited at most once.
Pheromone and heuristic information: The pheromone in the TSP problem represents the expectation of visiting city J directly after visiting city I. The heuristic information value is inversely proportional to the distance between city I and city J.
Solution construction: Each ant starts from a randomly selected city, and after each iteration, the ant adds a city that has not yet been visited to the solution. When all cities have been visited by ants, the construction of the solution terminates.
The ant colony algorithm has defects:
Ant colony algorithm can barely be used to solve small-scale TSP problems, and it takes a little time to find the optimal solution. However, if the problem is large, the performance of ant colony algorithm will be extremely low, or even stuck. So you can make improvements, like the Elite ant system.
The elite ant system is an improvement of the basic ant colony algorithm, which adds an enhancement method to the optimal path on the basis of the original AS pheromone updating principle. At the end of each pheromone update, the ant that has found the best path so far will add additional pheromones to the path. The introduction of such additional pheromone enhancement means in the elite ant system helps to better guide the bias of ant search and make the algorithm converge faster
Second, examples and codes
Content of Excel EXP12_3_2.xls:
function []=antvrptwwithjamminallcost()clear all; %清除所有变量 close all; %清图 clc ; %清屏 tic%开始计时。考虑拥堵时间段，以总费用作为比较目标%%SOLOMON问题数据
wenjian=textread('RXx201.txt');%%调用文件 changdu=size(wenjian,1);%%调用文件的长度 city_coordinate=zeros(changdu,2);%%需求点坐标 demands=zeros(1,changdu);%%需求 earlytime=zeros(1,changdu);%%需求点时间最早服务限制 windowtime=zeros(1,changdu);%%需求点时间限制 servicetime=zeros(1,changdu);%服务时间 for i=1:changdu city_coordinate(i,1)=(wenjian(i,2));%坐标 city_coordinate(i,2)=(wenjian(i,3));%坐标demands(i)=(wenjian(i,4));%需求量earlytime(i)=10;%(wenjian(i,5))-30;%时间限制,最早时间窗windowtime(i)=(wenjian(i,6))+60;%时间限制servicetime(i)=(wenjian(i,7));%服务时间 end
% city_coordinate %需求点坐标% demands%需求量% earlytime%%需求点时间最早服务限制% windowtime%时间窗限制% servicetime%服务时间
% city_coordinate=[42,59;6,17;37,19;22,76;28,11;21,16;12,65;35,18;38,29];%坐标 % demands=[0,90,40,60,70,70,40,20,40];%需求量% windowtime=[0,100,200,300,300,300,300,300,300];%时间窗限制% earlytime=[0,20,40,230,110,220,40,220,40];%服务时间% servicetime=[0,20,40,30,10,20,40,20,40];%服务时间
%%初始化变量m=40;% m 蚂蚁个数NC_max=10;% 最大允许运行次数Alpha=1;% Alpha 表征信息素重要程度的参数Beta=3;% Beta 表征启发式因子重要程度的参数Rho=0.25;% Rho 信息素蒸发系数Q=20;% Q 信息素增加强度系数vehiclecapacity=200;%车辆容量C=city_coordinate;n=size(C,1);%n需求点个数%%计算节点之间的距离D=zeros(n,n);%D距离矩阵 for i=1:n for j=1:n if i~=j D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5; else D(i,j)=eps; end D(j,i)=D(i,j); endend D;%距离矩阵 %%carspeed=60;%运输车辆正常行驶速度jamcarspeed=20;%拥堵时段的车辆行驶速度load_w=0;%车辆载重Eta=1./D;% %Eta为启发因子矩阵，这里设为距离的倒数Tau=ones(n,n);%%Tau为信息素矩阵Tabu=zeros(m,n);%存储并记录路径的生成，禁忌表NC=1;%当前运行次数G_best_route=[NC_max,n*2];%各代最佳路线G_best_length=inf.*ones(1,NC_max);%%各代最佳路线的长度length_ave=zeros(1,NC_max);%%各代路线的平均长度%iterbestlength=zeros(1,NC_max);%每次迭代最优距离的长度dividelength=2;%道路分阶段的长度%%定义计算碳排放的变量bestfuelconsumptionofcar=0;%该车辆的碳排放量；iterbestroute=zeros(1,n*2);%各代最佳路线dangqianroute=zeros(1,n*2);%当前的各代最佳路线bestroute=zeros(1,n*2);%最佳路线bestjuli=0;%最佳路线长度bijiaomubiao=0;%每次的比较目标minmubiao=0;qiyoujiage=7.5;%汽油价格fixfeeofvehicle=200;%汽车固定使用费用varityfeeofvehicle=2;%汽车变动使用费用：单位：元/单位距离carbonunitoffuel=2.32;%每升汽油的碳排放标准,单位：公斤unitcarbonfee=0.0528;%每公斤碳的收费价格renlifeeunit=0.5;%每分钟的人力资源成本besttravletime=0;%所有服务车辆的行驶时间bestcarfee=0;%所有服务车辆的费用，车辆固定费用+变动使用费用+人力成本bestcarbon=0;%%所有服务车辆的碳排放数量bestqiyoufee=0;%%所有服务车辆的汽油费用
Delta_Tau=zeros(n,n);%%信息素矩阵
%% %停止条件之一：达到最大迭代次数，停止while NC<=NC_max%循环控制 %如randi(2,3,[1 6])，就是产生一个2*3随机矩阵，这个矩阵的元素是区间[1 6]的随机数。 Tabu(:,1)=randi(1,m,1); factfuelconsumptionofcar=zeros(1,m);%每次只蚂蚁访问所有路径的所有油耗 totaltimeofvehicletravle=zeros(1,m);%每次只蚂蚁访问所有路径的所有时间 anttraveldistance=zeros(1,m);%每次只蚂蚁访问所有节点的距离 factantcarbon=zeros(1,m);%蚂蚁访问所有路径的视角产生的碳排放量 routejuli=0;%%路线距离 bestmubiao=0;%最优的比较目标%%m只蚂蚁的访问遍历情况 for i=1:m% disp('===========新的蚂蚁开始访问============')% disp('蚂蚁所用时间')% totaltimeofvehicletravle(i) visited=Tabu(i,:);%已经访问的需求点 visited=visited(visited>0); to_visit=setdiff(1:n,visited);%比较2个数组内不同的值，得出还没有访问的需求点 %%蚂蚁开始访问 load_w=0;%% vehicletime=0;%%车辆时间控制 dangqianvehiclevisittime=0; carcarbon=0;%蚂蚁访问所有路径的可能产生的碳排放量 fuelconsumptionofcar=0;%蚂蚁访问所有路径的可能产生的油耗量 j=1; while j<=n% disp('访问当前节点的具体时间，新的节点')% vehicletime if ~isempty(to_visit)%如果还有没有访问的需求点 to_visit ; %% 判断选择哪个需求点 for k=1:length(to_visit)%对每一个将要访问的需求点进行信息素计算 x(k)=(Tau(visited(end),to_visit(k))^Alpha)*(Eta(visited(end),to_visit(k))^Beta); %Tau为信息素矩阵，% Alpha 表征信息素重要程度的参数，%Eta为启发因子矩阵，这里设为距离的倒数，%% Beta 表征启发式因子重要程度的参数 end %%需求点选择方法 ww=rand;%随机生成一个概率 x=x/(sum(x)); xcum=cumsum(x); Select=find(xcum>=ww);%若计算的概率大于原来的就选择这条路线%要选择其中总概率大于等于某一个随机数，找到大于等于这个随机数的需求点的在未访问需求点中的位置 %%选择到的节点 % disp('访问当前节点')% to_visit(Select(1)) load_w=load_w+demands(to_visit(Select(1)));%车辆装载量计算 r=load_w/vehiclecapacity;%汽车载重与车辆容量的比例 Tabu(i,length(visited));%当前已经访问节点 to_visit(Select(1));%当前将要访问节点% disp('访问当前节点的距离')% D(Tabu(i,length(visited)),to_visit(Select(1)))%节点之间的距离 %开始该路段的时间计算 pathlength=D(Tabu(i,length(visited)),to_visit(Select(1))); %该路段长度 kk=ceil(pathlength/dividelength);%ceil表示向上取整，得出该路段分得的所有段数 currentcartime=zeros(1,kk);%%车辆行驶在路段K的当前行驶时间 %%拥堵时间判断 if kk==1 & dividelength>pathlength %第一个路段非常小 if (vehicletime>=60 & vehicletime<=120) || (vehicletime>=720 & vehicletime<=840) currentcartime(1)=(pathlength/jamcarspeed)*60;%拥堵时间段行驶的各个路段时间 currentfuel(1)=pathlength*(0.0617*jamcarspeed*jamcarspeed-7.8227*jamcarspeed+429.51)/1000;%车辆油耗 carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/jamcarspeed-1.33/jamcarspeed; else currentcartime(1)=(pathlength/jamcarspeed)*60;%正常行驶的各个路段时间 currentfuel(1)=pathlength*(0.0617*carspeed*carspeed-7.8227*carspeed+429.51)/1000;%车辆油耗 carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/carspeed-1.33/carspeed; end vehicletime= vehicletime+currentcartime(1); fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(1); carcarbon=carcarbon+currentfuel(1)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量 else currentcartime(1)=(dividelength/carspeed)*60;%正常行驶的各个路段时间 vehicletime= vehicletime+currentcartime(1); currentfuel(1)=dividelength*(0.0617*carspeed*carspeed-7.8227*carspeed+429.51)/1000;%车辆油耗 fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(1); carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/carspeed-1.33/carspeed; carcarbon=carcarbon+currentfuel(1)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量 end for ii=2:kk%计算每一段路段的行驶速度与行驶时间 if ii<kk if (vehicletime>=60 & vehicletime<=120) || (vehicletime>=720 & vehicletime<=840) currentcartime(ii)=dividelength/jamcarspeed*60;%正常时间段行驶的各个路段时间 currentfuel(ii)=dividelength*(0.0617*jamcarspeed*jamcarspeed-7.8227*jamcarspeed+429.51)/1000;%车辆油耗 carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/jamcarspeed-1.33/jamcarspeed; else currentcartime(ii)=dividelength/carspeed*60;%正常时间段行驶的各个路段时间 currentfuel(ii)=dividelength*(0.0617*carspeed*carspeed-7.8227*carspeed+429.51)/1000;%车辆油耗 carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/carspeed-1.33/carspeed; end vehicletime=vehicletime+currentcartime(ii); fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(ii); carcarbon=carcarbon+currentfuel(ii)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量 end if (ii==kk) %最后一个路段的各个路段时间 if (vehicletime>=60 & vehicletime<=120) || (vehicletime>=720 & vehicletime<=840) currentcartime(ii)=dividelength/jamcarspeed*60;%正常时间段行驶的各个路段时间 currentfuel(ii)=dividelength*(0.0617*jamcarspeed*jamcarspeed-7.8227*jamcarspeed+429.51)/1000;%车辆油耗 carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/jamcarspeed-1.33/jamcarspeed; else currentcartime(ii)=dividelength/carspeed*60;%正常时间段行驶的各个路段时间 currentfuel(ii)=dividelength*(0.0617*carspeed*carspeed-7.8227*carspeed+429.51)/1000;%车辆油耗 carbonxiuzhengbili=1.27+0.0614*r-0.0011*r*r-0.00235/carspeed-1.33/carspeed; end vehicletime=vehicletime+currentcartime(ii); fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(ii); carcarbon=carcarbon+currentfuel(ii)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量 end end% disp('+++++随机需求点：车辆行驶该路段的行驶时间')% vehicletime % vehicletime= vehicletime+servicetime(to_visit(Select(1))); %将要访问的需求点的服务时间计算 % disp('车辆行驶该路段的行驶时间++加上服务时间')% servicetime(to_visit(Select(1)))% vehicletime dangqianvehiclevisittime=dangqianvehiclevisittime+vehicletime; %%%容量和时间判断 if (load_w<=vehiclecapacity && dangqianvehiclevisittime<=windowtime(to_visit(Select(1))) ) && dangqianvehiclevisittime>=earlytime(to_visit(Select(1))) %if dangqianvehiclevisittime<=windowtime(to_visit(Select(1))) && dangqianvehiclevisittime>=earlytime(to_visit(Select(1))) Tabu(i,length(visited)+1)=to_visit(Select(1)); %访问的节点 totaltimeofvehicletravle(i)=totaltimeofvehicletravle(i)+vehicletime;%%蚂蚁i所有的车辆行驶时间% disp('+++++++++车辆的行驶时间之和+++++++++++++')% totaltimeofvehicletravle(i) anttraveldistance(i)=anttraveldistance(i)+pathlength;%%蚂蚁i所有车辆的行驶距离 factantcarbon(i)=factantcarbon(i)+carcarbon;%%蚂蚁i所有车辆的实际产生的碳排放 factfuelconsumptionofcar(i)=factfuelconsumptionofcar(i)+fuelconsumptionofcar; vehicletime=0;%服务时间归零 fuelconsumptionofcar=0; carcarbon=0; else% disp('该车辆服务完毕的的行驶时间')% vehicletime j=j-1;%没有选择到需求点，因此循环计数器减少一次 load_w=0;%车辆装载量归零 vehicletime=0;%服务时间归零 carcarbon=0;%碳排放归零 fuelconsumptionofcar=0; Tabu(i,length(visited)+1)=1; %返回了起点 dangqianvehiclevisittime=0; end end %%访问了需求点
% %if (load_w>vehiclecapacity) || ((dangqianvehiclevisittime>windowtime(to_visit(Select(1))))&&(dangqianvehiclevisittime<earlytime(to_visit(Select(1)))))% if (load_w>vehiclecapacity) || (dangqianvehiclevisittime>windowtime(to_visit(Select(1))))% % disp('该车辆服务完毕的的行驶时间')% % vehicletime% j=j-1;%没有选择到需求点，因此循环计数器减少一次% load_w=0;%车辆装载量归零% vehicletime=0;%服务时间归零% carcarbon=0;%碳排放归零% fuelconsumptionofcar=0;% Tabu(i,length(visited)+1)=1; %返回了起点 % dangqianvehiclevisittime=0;% else % % Tabu(i,length(visited)+1)=to_visit(Select(1)); %访问的节点% totaltimeofvehicletravle(i)=totaltimeofvehicletravle(i)+vehicletime;%%蚂蚁i所有的车辆行驶时间% % disp('+++++++++车辆的行驶时间之和+++++++++++++')% % totaltimeofvehicletravle(i)% anttraveldistance(i)=anttraveldistance(i)+pathlength;%%蚂蚁i所有车辆的行驶距离% factantcarbon(i)=factantcarbon(i)+carcarbon;%%蚂蚁i所有车辆的实际产生的碳排放% factfuelconsumptionofcar(i)=factfuelconsumptionofcar(i)+fuelconsumptionofcar;% vehicletime=0;%服务时间归零% fuelconsumptionofcar=0;% carcarbon=0;% end
% disp('实际车辆行驶的时间')% vehicletime% totaltimeofvehicletravle(i) j=j+1; visited=Tabu(i,:);%已经访问的需求点 visited=visited(visited>0); to_visit=setdiff(1:n,visited);%比较2个数组内不同的值，得出还没有访问的需求点 x=[];%清空没有访问的需求点的信息 if visited(end)~=1 Tabu(i,1:(length(visited)+1))=[visited,1]; end end % disp('车辆实际行驶的总时间')% totaltimeofvehicletravle(i)% disp('车辆实际行驶的路线')% Tabu(i,:) end %%全部M只蚂蚁循环访问控制
factfuelconsumptionofcar;factantcarbon;totaltimeofvehicletravle;anttraveldistance;Tabu;
%%计算m只蚂蚁的各自的耗油量、行驶时间、行驶距离带来的总费用vehiclenumofeachant=zeros(1,m);%各只蚂蚁的总费用车辆数量anttotalcost=zeros(1,m);%各只蚂蚁的总费用初始化for i=1:m nn=Tabu(i,:); for j=2:length(nn) if nn(j)==1 vehiclenumofeachant(i)=vehiclenumofeachant(i)+1; end end anttotalcost(i)=anttraveldistance(i)*varityfeeofvehicle+vehiclenumofeachant(i)*fixfeeofvehicle+totaltimeofvehicletravle(i)*renlifeeunit+factfuelconsumptionofcar(i)*qiyoujiage+factantcarbon(i)*unitcarbonfee;%endbestcost=min(anttotalcost);%%找出总费用最优的路线weizhi=find(bestcost==anttotalcost); %%找出总费用最优的路线的位置kk=Tabu(weizhi,:);for x=1:length(kk)iterbestroute(x)=kk(x);end
%%每次最好的计算结果记录if NC>=2 for x=1:length(kk) Tabu(1,x)=iterbestroute(x); end enddangqianbestroute=iterbestroute(iterbestroute>0);%length_ave(NC)=mean(anttotalcost);%%以总费用作为目标的平均值length_ave(NC)=bestcost;%%以总费用作为目标的平均值
bijiaomubiao=bestcost; %%本次程序的比较目标if NC==1 bestfuelconsumptionofcar=factfuelconsumptionofcar(weizhi); bestroute=dangqianbestroute; bestjuli=anttraveldistance(weizhi); minmubiao=bijiaomubiao; besttravletime=totaltimeofvehicletravle(weizhi); bestcarbon=factantcarbon(weizhi); else if bijiaomubiao<minmubiao bestfuelconsumptionofcar=factfuelconsumptionofcar(weizhi); bestroute=dangqianbestroute; bestjuli=anttraveldistance(weizhi); minmubiao=bijiaomubiao; besttravletime=totaltimeofvehicletravle(weizhi); bestcarbon=factantcarbon(weizhi); endend
%% 第四步记录本代各种参数L=zeros(m,1);for i=1:m MM=Tabu(i,:); R=MM(MM>0); for j=1:(length(R)-1) L(i)=L(i)+D(R(j),R(j+1)); end end NC=NC+1 %迭代计数器%% 全局信息素更新%Delta_Tau=zeros(n,n); for i=1:m MM=Tabu(i,:); R=MM(MM>0); cd=length(R); for j=1:cd Tabu(i,j)=R(j); end for j=1:(cd-1) %Delta_Tau(R(j),R(j+1))=Delta_Tau(R(j),R(j+1))+Q/L(i); Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i); %此次循环在路径（i，j）上的信息素增量 end %此次循环在整个路径上的信息素增量 end Tau=(1-Rho).*Tau+Delta_Tau; %考虑信息素挥发，更新后的信息素
%% ????????? Tabu=zeros(m,n); load_w=0;end
%%输出最优路线的相关信息 minmubiao bestfuelconsumptionofcar bestroute bestjuli besttravletime bestcarbon
%% 画图画出路线subplot(1,2,1) plot([C(bestroute,1)],[C(bestroute,2)],'-*')%最优路径 subplot(1,2,2) %绘制第二个子图形% plot(iterbestlength)%各代的最短距离% hold onplot(length_ave,'r')%各代的平均距离% title('平均距离和最短距离') %标题tocend
Copy the code
The complete code to download www.cnblogs.com/ttmatlab/p/…