function [G_k,G_best_route,G_best_length,K,best_route,best_length,length_ave]=VRP(C,NC_max,m,Alpha,Beta,Rho,Q,W,t) %#ok<INUSD,STOUT>
%%
%--------------------------------------------------------------------------
%%参数说明
%%G_K:各代最优车辆数目
%%G_best_route:各代最优路线
%%G_best_length:各代最优路线长度
%%K:最优车辆数目
%%best_route:最优路线
%%best_length:最优路线长度
%%length_ave:各代平均长度
%%C:DC和工厂的坐标
%%NC_max:最大迭代次数
%%m:蚂蚁数目
%%Alpha:重要度系数
%%Beta:能见度系数
%%Rho:挥发度系数
%%Q:信息更新参数
%%W:车辆载重量
%%
%%第一步 初始化变量和系数 m=60;Alpha=1;Beta=1;gama=2;Rho=0.15;NC_max=50;Q=15;W=9;qq=0.05; %#ok<NASGU>
C=[0 0
0 -1
0 3
-2 -2
-3 -3
3 -1
-4 0
-4 -1
1 -2
1 -1
1 3
3 4
-3 0
2 0
1 -3
2 -1
2 1
1 -4
-3 2
-1 -1]; %#ok<NASGU>
t=[0 1.5 1.8 2 0.8 1.5 1.0 2.5 3.0 1.7 0.6 0.2 2.4 1.9 2.0 0.7 0.5 2.2 3.1 0.1]; %构造仓库和工厂等之间的距离矩阵 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);
end
end %构造节省量矩阵
U=zeros(n,n);%U表示工厂之间的连接和其与仓库之间连接能节省的距离
for i=1:n
for j=1:n
if i~=j
U(i,j)=D(i,1)+D(j,1)-D(i,j);
else
U(i,j)=eps;
end
U(j,i)=U(i,j);
end
end
load_w=0;
Eta=1./D;%Eta为启发因子,这里设为距离的倒数
Tau=ones(n,n);%Tau为信息素矩阵
Tabu=zeros(m,n+20);%存储并记录路径的生成
NC=1;%迭代计数器
G_best_route=[NC_max,n+20];%各代最佳路线
G_best_length=inf.*ones(NC_max,1);%各代最佳路线的长度
length_ave=zeros(NC_max,1);%各代路线的平均长度 %%第二步,把蚂蚁放到DC内
while NC<=NC_max%停止条件之一:达到最大迭代次数
Tabu(:,1)=randint(m,1,[1,1]);
%%第三步,m只蚂蚁按照要求的方法选择工厂,并完成周游
for i=1:m
visited=Tabu(i,:);
visited=visited(visited>0);
to_visit=setdiff(1:n,visited);
c_temp=length(to_visit);
j=1;
while j<=n
if ~isempty(to_visit)
%visit_thisant=[];
%while to_visit~=[1]
%for ee=1:(2*c_temp)
%while c_temp>=2
%if length(to_visit)>1
%visited=Tabu(i,:);
%to_visit=setdiff(1:n,visited);
%to_visit=[1,to_visit];
%c_temp=length(to_visit);
%visit_thisant=[];
% if c_temp~=1
%for b=1:c_temp
%if (load_w+t(to_visit(b)))<=W
%visit_thisant=[visit_thisant,to_visit(b)];
%end
%end
%jj=visited(visited>0);
%% 按照规则选下一个工厂或者是回到仓库
%for k=1:length(visit_thisant)
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);