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

Posted 2021-06-26 博主QQ2449341593

一、原理简介

蚁群算法的基本原理

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

完整代码或者仿真咨询QQ1575304183