可以直接运行求解VRP的matlab程序，

function []=antvrptw（）
clc
clear
tic%开始计时
city_coordinate=[33042312;36391315;41772244;37121399;34881535;33261556;32381229;41961044;4312790;4386170;
                 30071970;25621756;27881491;23811676;1332695;37151678;39182179;40612370;37802212;36762578;
                 30291838];%坐标           
nodedemand=[09040607070402040706040404020805030505050];%需求量
timewindow=[0290240260270270240220240270260240340340320380350230250250350];%时间窗限制
servicewindow=[09040607070402040706040404020805030505050];%服务时间

m=30;% m 蚂蚁个数
Alpha=3;% Alpha 表征信息素重要程度的参数
Beta=2;% Beta 表征启发式因子重要程度的参数
gama=2;
Rho=0.2;% Rho 信息素蒸发系数
NC_max=50;
Q=2;% Q 信息素增加强度系数
W=300;
qq=0.1;
load_w=0;
s=0;%??????
C=city_coordinate;
demands=nodedemand;
windowtime=timewindow;
servicetime=servicewindow;
vehicletime=0;

n=size（C1）;%n????????????? 
D=zeros（nn）;%D???????????? 
for i=1:n 
   for j=1:n 
     if i~=j 
         D（ij）=（（C（i1）-C（j1））^2+（C（i2）-C（j2））^2）^0.5; 
     else 
        D（ij）=eps; 
     end 
    D（ji）=D（ij）; 
   end 
end 
Eta=1./D;%Eta??????????????? 
Tau=ones（nn）;%Tau?????? 
Tabu=zeros（mn+32）;%?????????? 
NC=1;%????? 
G_best_route=[NC_maxn+32];%?????? 
G_best_length=inf.*ones（NC_max1）;%????????? 
length_ave=zeros（NC_max1）;%?????????

%%?????????DC?
while NC<=NC_max%??????????????? 
    %Tabu（:1）=randint（m1[11]）;%如randint（23[1 6]），就是产生一个2*3随机矩阵，这个矩阵的元素是区间[1 6]的随机数。
    Tabu（:1）=randi（1m1）;
%%????m???????????????????? 
for i=1:m
    visited=Tabu（i:）;
    visited=visited（visited>0）;
    to_visit=setdiff（1:nvisited）;
    c_temp=length（to_visit）;
    j=1;
       while j<=n
          if ~isempty（to_visit）
   %% ?????????????????
    for k=1:length（to_visit）
        x（k）=（Tau（visited（end）to_visit（k））^Alpha）*（Eta（visited（end）to_visit（k））^Beta）;%*（U（visited（end）to_visit（k））^gama）;
          % x（k）=（Tau（visited（end）to_visit（k））^Alpha）*（Eta（visited（end）to_visit（k））^Beta）*（U（visited（end）to_visit（k））^gama）;
    end   
         ww=rand;
    if ww         Select=find（max（x））;
    el