一、原理简介
蚁群算法的基本原理
1、蚂蚁在路径上释放信息素。
2、碰到还没走过的路口，就随机挑选一条路走。同时，释放与路径长度有关的信息素。
3、信息素浓度与路径长度成反比。后来的蚂蚁再次碰到该路口时，就选择信息素浓度较高路径。
4、最优路径上的信息素浓度越来越大。
5、最终蚁群找到最优寻食路径。
蚁群走过较短的那一侧的蚂蚁数数量会多于较长那一侧的，所以留下的信息素就会多，渐渐的蚂蚁就只走较短的那一侧了。
蚁群算法对TSP的求解主要有两大步骤：（TSP问题就是要找到最短哈密尔顿回路）
1、路径构建
AS中的随机比例规则；对于每只蚂蚁k，路径记忆向量R^k.按照访问顺序记录了所有k已经经过的城市序号。设蚂蚁k当前所在城市为i,则其选择城市j作为下一个访问对象的概率为:
2、信息素更新
这里m是蚂蚁个数, ρ是信息素的蒸发率，规定0≤ ρ≤1，在AS中通常设置为 ρ =0.5，Δτij是第k只蚂蚁在它经过的边上释放的信息素量，它等于蚂蚁k本轮构建路径长度的倒数。C^k表示路径长度，它是R^k中所有边的长度和。
构建图：构建图与问题描述图是一致的，成份的集合C对应着点的集合（即：C=N），连接对应着边的集合（即L=A）,且每一条边都带有一个权值，代表点i和j之间的距离。
约束条件：所有城市都要被访问且每个城市最多只能被访问一次。
信息素和启发式信息：TSP 问题中的信息素表示在访问城市i后直接访问城市j的期望度。启发式信息值一般与城市i和城市j的距离成反比。
解的构建：每只蚂蚁最初都从随机选择出来的城市出发，每经过一次迭代蚂蚁就向解中添加一个还没有访问过的城市。当所有城市都被蚂蚁访问过之后，解的构建就终止。
蚁群算法存在缺陷：
蚁群算法在解决小规模TSP问题是勉强能用，稍加时间就能发现最优解，但是若问题规模很大，蚁群算法的性能会极低，甚至卡死。所以可以进行改进，例如精英蚂蚁系统。
精英蚂蚁系统是对基础蚁群算法的一次改进，它在原AS信息素更新原则的基础上增加了一个对至今最优路径的强化手段。在每轮信息素更新完毕后，搜索到至今最优路径的那只蚂蚁将会为这条路径添加额外的信息素。精英蚂蚁系统引入 这种额外的信息素强化手段有助于更好的引导蚂蚁搜索的偏向，使算法更快收敛
二、实例及代码
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:changducity_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:nfor j=1:nif i~=jD(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;elseD(i,j)=eps;endD(j,i)=D(i,j);endendD;%距离矩阵%%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('访问当前节点的具体时间，新的节点')% vehicletimeif ~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;elsecurrentcartime(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;endvehicletime= vehicletime+currentcartime(1);fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(1);carcarbon=carcarbon+currentfuel(1)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量elsecurrentcartime(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;%%这只蚂蚁的碳排放量endfor ii=2:kk%计算每一段路段的行驶速度与行驶时间if ii<kkif (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;elsecurrentcartime(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;endvehicletime=vehicletime+currentcartime(ii);fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(ii);carcarbon=carcarbon+currentfuel(ii)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量endif (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;elsecurrentcartime(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;endvehicletime=vehicletime+currentcartime(ii);fuelconsumptionofcar=fuelconsumptionofcar+currentfuel(ii);carcarbon=carcarbon+currentfuel(ii)*carbonxiuzhengbili*carbonunitoffuel;%%这只蚂蚁的碳排放量endend% disp('+++++随机需求点：车辆行驶该路段的行驶时间')% vehicletime%vehicletime= vehicletime+servicetime(to_visit(Select(1))); %将要访问的需求点的服务时间计算% disp('车辆行驶该路段的行驶时间++加上服务时间')% servicetime(to_visit(Select(1)))% vehicletimedangqianvehiclevisittime=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('该车辆服务完毕的的行驶时间')% vehicletimej=j-1;%没有选择到需求点，因此循环计数器减少一次load_w=0;%车辆装载量归零vehicletime=0;%服务时间归零carcarbon=0;%碳排放归零fuelconsumptionofcar=0;Tabu(i,length(visited)+1)=1; %返回了起点dangqianvehiclevisittime=0;endend %%访问了需求点% %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)~=1Tabu(i,1:(length(visited)+1))=[visited,1];endend% 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:mnn=Tabu(i,:);for j=2:length(nn)if nn(j)==1vehiclenumofeachant(i)=vehiclenumofeachant(i)+1;endendanttotalcost(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>=2for x=1:length(kk)Tabu(1,x)=iterbestroute(x);endenddangqianbestroute=iterbestroute(iterbestroute>0);%length_ave(NC)=mean(anttotalcost);%%以总费用作为目标的平均值length_ave(NC)=bestcost;%%以总费用作为目标的平均值bijiaomubiao=bestcost; %%本次程序的比较目标if NC==1bestfuelconsumptionofcar=factfuelconsumptionofcar(weizhi);bestroute=dangqianbestroute;bestjuli=anttraveldistance(weizhi);minmubiao=bijiaomubiao;besttravletime=totaltimeofvehicletravle(weizhi);bestcarbon=factantcarbon(weizhi);else if bijiaomubiao<minmubiaobestfuelconsumptionofcar=factfuelconsumptionofcar(weizhi);bestroute=dangqianbestroute;bestjuli=anttraveldistance(weizhi);minmubiao=bijiaomubiao;besttravletime=totaltimeofvehicletravle(weizhi);bestcarbon=factantcarbon(weizhi);endend%% 第四步记录本代各种参数L=zeros(m,1);for i=1:mMM=Tabu(i,:);R=MM(MM>0);for j=1:(length(R)-1)L(i)=L(i)+D(R(j),R(j+1));endendNC=NC+1 %迭代计数器%% 全局信息素更新%Delta_Tau=zeros(n,n);for i=1:mMM=Tabu(i,:);R=MM(MM>0);cd=length(R);for j=1:cdTabu(i,j)=R(j);endfor 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%此次循环在整个路径上的信息素增量endTau=(1-Rho).*Tau+Delta_Tau; %考虑信息素挥发，更新后的信息素%% ?????????Tabu=zeros(m,n);load_w=0;end%%输出最优路线的相关信息minmubiaobestfuelconsumptionofcarbestroutebestjulibesttravletimebestcarbon%% 画图画出路线subplot(1,2,1)plot([C(bestroute,1)],[C(bestroute,2)],'-*')%最优路径subplot(1,2,2) %绘制第二个子图形% plot(iterbestlength)%各代的最短距离% hold onplot(length_ave,'r')%各代的平均距离% title('平均距离和最短距离') %标题tocend
