【VRP问题】基于蚁群算法求解带时间窗车辆调度问题

一、原理简介

蚁群算法的基本原理

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