一、原理简介
蚁群算法的基本原理
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: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
