基于模拟退火和蚁群算法求解有载重限制的车辆路径优化问题
VRP问题介绍
车辆路径优化问题（Vehicle RoutingProblem，VRP）由Dantzig于1959年提出，是一个经典的NP 问题。问题描述：车辆从配送中心出发，依次对随机分布的多个配送点执行配送任务，每个配送点有且仅有一辆车辆通过，车辆完成配送任务后返回配送中心。要求在符合约束条件的情况下，找到一条最佳路径，代价最低。根据具体约束条件的不同，VRP问题可以分为不同的类型，如受载重量约束的CVRP问题(Capacitated Vehicle Routingroblem,CVRP)等。
可以选择精确算法和启发式算法求解VRP问题。精确算法能够求得VRP问题的最优解，但随着配送点的增加，问题复杂度呈指数级增长，在大规模配送点的情况下并不适用。
目前启发式算法已成为解决该类问题的主要方法，原理是通过计算机程序仿真自然界中的群体智慧来解决实际问题。典型算法有蚁群算法(Ant ColonyOptimization，ACO)、遗传算法(Genetic Algorithms，GA)、禁忌搜索算法(Tabu Search，TS)和神经网络算法(Neural Network，NN)等。
本文针对蚁群算法加速收敛和停滞现象的矛盾，提出一种基于模拟退火机制的蚁群算法。在高温阶段以高的概率接收候选集内较差的解加入更新集,随着温度的降低，路径上的信息素分布趋于集中，可以使算法快速收敛，在温度控制机制上采用回火策略。
CVRP数学模型
设配送中心有k辆车，向n个配送点配送货物，配送模型用一个加权图G(V,E)来表示，其中V(v0,v1,…,vn)表示配送点集，v0代表配送中心，其余为客户点，E={(vi,vj)|vi,vj∈V;i≠j}代表从点i到点j的行驶路径集合。设每辆车的载重量为D，各客户点的需求为qi(i∈V)。以c表示车辆行驶距离的单位成本，dij表示从配送点i到配送点j的行驶距离，车辆数为m，可建立CVRP数学模型如下：
模拟退火蚁群混合思路：
设系统温度为Tc,初始温度Tc(0)=T0，第t次迭代时的温度Tc(t)=Tt.当m只蚂蚁完成一轮搜索时，取得解集称为候选集{bk(sk,lk)1≤k≤m}，式中：bk(sk,lk)是蚂蚁k取得的解；sk是经过的城市路径；lk是路径长度.候选集中最优个体记为bl-min(sl-min,ll-min)，当前算法得到的最优个体记为bmin(smin,lmin)，其中：sl-min为局部最优解中蚂蚁走过的路径；ll-min为局部最优解中最短路径长度。
采用模拟退火机制由候选集生成更新集，更新集记作{b′k(sk，lk)1≤k≤m}，逐一计算候选集内的各个蚂蚁的解，根据模拟退火机制决定候选集内的解bk(sk,lk)是否加入更新集，即：
若pt=1，则接收bk(sk,lk)作为更新集内的b′k(sk，lk)；否则产生一个[0,1)区间的均匀分布的随数ξ，若pt>ξ，则bk进入更新集，否则bl-min进入更新集。对所有候选集内的个体进行筛选形成更新集{b′k(sk，lk)1≤k≤m}。
温度控制策略与回火机制
在完成一轮搜索与路径信息素更新后，进行降温，降温函数为：
式中a为降温系数。可以看出，当系统温度Tc高时，以高的概率接收较差解加入更新集，使路径上的信息素分布更加广泛，避免出现早熟；随着Tc的降低，接受概率减少，对较差解接收的概率也变小，使路径上的信息素分布集中，算法得以快速收敛。在本文算法中，鼓励过早收敛，采用较小的a，这样可以缩短退火过程找到更好解的时间。当然，这会增大落入局部极值点的风险。为了解决这个问题，算法中引入回火机制，回火的温度范围为[Tmax,Tmin]，当Tc<Tmin时，将系统温度升高至Tmax，重新进行算法迭代；当再次Tc<Tmin时，这一过程称为回火过程。可以根据设定采取多次回火，回火次数记为H，最大回火次数记为Hmax。
参考文献
[1]刘波,蒙培生.采用基于模拟退火的蚁群算法求解旅行商问题[J].华中科技大学学报(自然科学版).
[2]刘希洋,赵建民,徐慧英,朱信忠.基于改进型蚁群算法求解车辆路径优化问题的研究[J].计算机时代.
部分程序
% Author: 怡宝2号
% Use: 用模拟退火蚁群算法求解有容量限制的车辆路径问题；
% 坐标视自己的具体情况而定。
% Illustrate:输入变量（must）：
% coordinate：要求的相应点的坐标，第一行是点的横坐标，第二行是点的纵坐标
% demand（1*n），n表示工厂的个数；t表示各个工厂要求运送的货物的重量
% W:每辆车限制的载重量
% Can—changedparameter： m：蚁群的规模
% NC_max：蚁群的最大迭代代数
% tmax：退火初始温度
% tmin：退火最低温度
% down：退火算法的降火系数
%
% 输出： best_route：最短路程对应的路径
% best_length ：最短路程
% attention： 如有疑问请咨询QQ：778961303
closeall
clearall
clc
%------参数初始化-----%
[m,Alpha,Beta,Gama,Rho,NC_max,Q,W,Prab]=initial();
num_client= 22; %客户的个数
%---------------模拟退火参数-----
tmax=20; %初始温度
tmin=1; %最低温度
down=0.85; %降温系数
t =tmax; %模拟退火温度初始化
coordinate=[
0 0 0 -2 -3 3 -4 -4 1 1 1 3 -3 2 1 2 2 1 -3 -1 4 4;
0 -1 3 -2 -3 -1 0 -1 -2 -1 3 4 0 0 -3 -1 1 -4 2 -1 2 -2]; %各个工厂或者客户的坐标（地址），第一行是x坐标，第二行是y坐标
coordinate=coordinate';
demand=[01.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 2 1.1]; %各个客户点的需求
demand=demand';
n=size(coordinate,1); %n表示问题的规模（城市个数）；取C的第一个值=n
Distance=zeros(n,n); %D表示各个城市之间的距离
g_route=[NC_max,n+20]; %各代最佳路线
temp_best_length=inf.*ones(NC_max,1); %各代最佳路线的长度
length_ave=zeros(NC_max,1); %各代路线的平均长度
temp_pra=0;
%---------------绘制仓库和客户点--------------------%
figure(1)
plot(coordinate(1,1),coordinate(1,2),'r-*');
text(coordinate(1,1),coordinate(1,2),[ ' ''仓库'] )
holdon
fori=2:n
plot(coordinate(i,1),coordinate(i,2),'-*');
text( coordinate(i,1),coordinate(i,2),[' ' num2str(i-1)] )
hold on
end
%--------------------求距离--------------------
fori=1:n
for j=1:n
if i~=j
Distance(i,j)=((coordinate(i,1)-coordinate(j,1))^2+(coordinate(i,2)-coordinate(j,2))^2)^0.5;
else
Distance(i,j)=eps;
end
Distance(j,i)=Distance(i,j);
end
end
% 构造工厂之间比分别从工厂出发时所能节省的路程
U=zeros(n,n);
fori=1:n
for j=1:n
if i~=j
U(i,j)=Distance(i,1)+Distance(j,1)-Distance(i,j);
else
U(i,j)=eps;
end
U(j,i)=U(i,j);
end
end
load=0; %
Eta=1./Distance; %Eta为启发因子，离散模型设为距离的倒数
Tau=ones(n,n); %Tau为信息素矩阵，行：前一个城市的标号；列：后一个城市的标号。
Tabu=zeros(m,n+15); %存储并记录路径的生成
NC=1; %迭代计数器
%------------------迭代开始------------------%
whileNC<=NC_max
Tabu(:,1)=1; %每只蚂蚁都是从仓库（及节点1）开始出发
for i=1:m
visited=Tabu(i,:);
visited=visited(visited>0);
to_visit=setdiff(1:n,visited); %在集合论中，to_visit = (1:n) - visited
c_temp=length(to_visit);
j=1;
%------------------------------------------------------------%
% 第一步按概率选择城市，第二步根据车子的载重量限制判定是否选择这个城市%
%------------------------------------------------------------%
while j<=n %参观了的城市数目
if ~isempty(to_visit) %判断to_visit不是空，是否还有城市没有访问到
%-----------------求解visited的最后一个城市，到所有还没有参观的城市的概率-----------------%
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);
end
temp_pra=rand;
if temp_pra<Prab %如果满足伪随机选择概率，直接选择概率最大的，这样收敛更快
Select=find(max(x));
else
x=x/(sum(x)); %用轮盘赌法选择下一个城市
xcum=cumsum(x);
Select=find(xcum>=rand);
end
if isempty(Select) %如果选不了下一个城市了，就回到1点
Select=1;
load=load+demand(Select); %select=1；t（select）=0
else
load=load+demand(to_visit(Select(1))); %select表示城市标号，to_visit(select)才表示城市
end
end
visited=Tabu(i,:);
visited=visited(visited>0); %已经拜访过的城市
to_visit=setdiff(1:n,visited);
x=[];
end
load=0;
end
%------------------求蚁群m只蚂蚁路径的路程------------------%
L=zeros(m,1); %共有m只蚂蚁
for i=1:m
temp=Tabu(i,:);
R=temp(temp>0);
for j=1:(length(R)-1)
L(i)=L(i)+Distance(R(j),R(j+1));
end
end
minL = min(L);
%%------------------引入模拟退火，产生更新集-------------------%
[Tabu,L]=SA(Tabu,L,m,n,demand,W,Distance,minL,t);
%---------模拟退火引入完成------------%
temp_best_length(NC)=min(L);
pos=find(L==temp_best_length(NC));
g_route(NC,1:length(Tabu(pos(1),:)))=Tabu(pos(1),:);
%% 画图，做成视频
figure(2)
DrawRoute(coordinate,temp);
title(['进化第',num2str(NC),'代']);
PICTURE(NC)=getframe;
hold off
%%2-opt优化完毕
length_ave(NC)=mean(L);
NC=NC+1 %显示进化的代数
%------------------更新信息素-------------------%
Delta_Tau=zeros(n,n);
for i=1:m
temp=Tabu(i,:);
R=temp(temp>0);
for j=1:(length(R)-1)
Delta_Tau(R(j),R(j+1))=Delta_Tau(R(j),R(j+1))+Q/L(i);
end
end
Tau=(1-Rho).*Tau+Delta_Tau;
%----------------更新模拟退后的温度---------------
t=t*down;
if t<tmin
t = tmax;
end
%------------------路径矩阵清零-------------------%
Tabu=zeros(m,n);
load=0;
end
%-------------------输出结果----------------------%
Pos=find(temp_best_length==min(temp_best_length)); %找到最短路径对应的行
best_route=g_route(Pos(1),:); %最短的路径
best_length=temp_best_length(Pos(1)); %最短路径的长度
%画图
best_route=best_route(best_route>0);
disp(['最短路线为：',num2str(best_route)])
disp(['最短路线对应的线路长度为：',num2str(best_length)])
%-------------------绘制迭代优化结果图----------------------%
%figure()
% p= plot(temp_best_length);
%set(p,'Color','red','LineWidth',2);
%xlabel('迭代代数')
%ylabel('目标函数值')
%title('进化过程')
%-------------------绘制迅游过程图----------------------%
figure(3)
plot(coordinate(1,1),coordinate(1,2),'r-*');
text( coordinate(1,1),coordinate(1,2),[' ''仓库'] )
hold on
fori=2:n
plot(coordinate(i,1),coordinate(i,2),'-*');
text( coordinate(i,1),coordinate(i,2),[' ' num2str(i-1)] )
hold on
end
figure(3)
%plot([coordinate(best_route,1)],[coordinate(best_route,2)],'r-*')
DrawRoute(coordinate,best_route);
holdon
PICTURE(NC_max)= getframe;
%制作视频
movie2avi(PICTURE,'SA-ACO-VRP.avi','fps',3,'quality',90)%输出视频
%movie2avi(PICTURE,'SA-ACO.avi', 'compression', 'None')%输出视频
结果：如有疑问，请咨询qq：778961303。
VRP问题介绍
车辆路径优化问题（Vehicle RoutingProblem，VRP）由Dantzig于1959年提出，是一个经典的NP 问题。问题描述：车辆从配送中心出发，依次对随机分布的多个配送点执行配送任务，每个配送点有且仅有一辆车辆通过，车辆完成配送任务后返回配送中心。要求在符合约束条件的情况下，找到一条最佳路径，代价最低。根据具体约束条件的不同，VRP问题可以分为不同的类型，如受载重量约束的CVRP问题(Capacitated Vehicle Routingroblem,CVRP)等。
可以选择精确算法和启发式算法求解VRP问题。精确算法能够求得VRP问题的最优解，但随着配送点的增加，问题复杂度呈指数级增长，在大规模配送点的情况下并不适用。
目前启发式算法已成为解决该类问题的主要方法，原理是通过计算机程序仿真自然界中的群体智慧来解决实际问题。典型算法有蚁群算法(Ant ColonyOptimization，ACO)、遗传算法(Genetic Algorithms，GA)、禁忌搜索算法(Tabu Search，TS)和神经网络算法(Neural Network，NN)等。
本文针对蚁群算法加速收敛和停滞现象的矛盾，提出一种基于模拟退火机制的蚁群算法。在高温阶段以高的概率接收候选集内较差的解加入更新集,随着温度的降低，路径上的信息素分布趋于集中，可以使算法快速收敛，在温度控制机制上采用回火策略。
CVRP数学模型
设配送中心有k辆车，向n个配送点配送货物，配送模型用一个加权图G(V,E)来表示，其中V(v0,v1,…,vn)表示配送点集，v0代表配送中心，其余为客户点，E={(vi,vj)|vi,vj∈V;i≠j}代表从点i到点j的行驶路径集合。设每辆车的载重量为D，各客户点的需求为qi(i∈V)。以c表示车辆行驶距离的单位成本，dij表示从配送点i到配送点j的行驶距离，车辆数为m，可建立CVRP数学模型如下：
模拟退火蚁群混合思路：
设系统温度为Tc,初始温度Tc(0)=T0，第t次迭代时的温度Tc(t)=Tt.当m只蚂蚁完成一轮搜索时，取得解集称为候选集{bk(sk,lk)1≤k≤m}，式中：bk(sk,lk)是蚂蚁k取得的解；sk是经过的城市路径；lk是路径长度.候选集中最优个体记为bl-min(sl-min,ll-min)，当前算法得到的最优个体记为bmin(smin,lmin)，其中：sl-min为局部最优解中蚂蚁走过的路径；ll-min为局部最优解中最短路径长度。
采用模拟退火机制由候选集生成更新集，更新集记作{b′k(sk，lk)1≤k≤m}，逐一计算候选集内的各个蚂蚁的解，根据模拟退火机制决定候选集内的解bk(sk,lk)是否加入更新集，即：
若pt=1，则接收bk(sk,lk)作为更新集内的b′k(sk，lk)；否则产生一个[0,1)区间的均匀分布的随数ξ，若pt>ξ，则bk进入更新集，否则bl-min进入更新集。对所有候选集内的个体进行筛选形成更新集{b′k(sk，lk)1≤k≤m}。
温度控制策略与回火机制
在完成一轮搜索与路径信息素更新后，进行降温，降温函数为：
式中a为降温系数。可以看出，当系统温度Tc高时，以高的概率接收较差解加入更新集，使路径上的信息素分布更加广泛，避免出现早熟；随着Tc的降低，接受概率减少，对较差解接收的概率也变小，使路径上的信息素分布集中，算法得以快速收敛。在本文算法中，鼓励过早收敛，采用较小的a，这样可以缩短退火过程找到更好解的时间。当然，这会增大落入局部极值点的风险。为了解决这个问题，算法中引入回火机制，回火的温度范围为[Tmax,Tmin]，当Tc<Tmin时，将系统温度升高至Tmax，重新进行算法迭代；当再次Tc<Tmin时，这一过程称为回火过程。可以根据设定采取多次回火，回火次数记为H，最大回火次数记为Hmax。
参考文献
[1]刘波,蒙培生.采用基于模拟退火的蚁群算法求解旅行商问题[J].华中科技大学学报(自然科学版).
[2]刘希洋,赵建民,徐慧英,朱信忠.基于改进型蚁群算法求解车辆路径优化问题的研究[J].计算机时代.
部分程序
% Author: 怡宝2号
% Use: 用模拟退火蚁群算法求解有容量限制的车辆路径问题；
% 坐标视自己的具体情况而定。
% Illustrate:输入变量（must）：
% coordinate：要求的相应点的坐标，第一行是点的横坐标，第二行是点的纵坐标
% demand（1*n），n表示工厂的个数；t表示各个工厂要求运送的货物的重量
% W:每辆车限制的载重量
% Can—changedparameter： m：蚁群的规模
% NC_max：蚁群的最大迭代代数
% tmax：退火初始温度
% tmin：退火最低温度
% down：退火算法的降火系数
%
% 输出： best_route：最短路程对应的路径
% best_length ：最短路程
% attention： 如有疑问请咨询QQ：778961303
closeall
clearall
clc
%------参数初始化-----%
[m,Alpha,Beta,Gama,Rho,NC_max,Q,W,Prab]=initial();
num_client= 22; %客户的个数
%---------------模拟退火参数-----
tmax=20; %初始温度
tmin=1; %最低温度
down=0.85; %降温系数
t =tmax; %模拟退火温度初始化
coordinate=[
0 0 0 -2 -3 3 -4 -4 1 1 1 3 -3 2 1 2 2 1 -3 -1 4 4;
0 -1 3 -2 -3 -1 0 -1 -2 -1 3 4 0 0 -3 -1 1 -4 2 -1 2 -2]; %各个工厂或者客户的坐标（地址），第一行是x坐标，第二行是y坐标
coordinate=coordinate';
demand=[01.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 2 1.1]; %各个客户点的需求
demand=demand';
n=size(coordinate,1); %n表示问题的规模（城市个数）；取C的第一个值=n
Distance=zeros(n,n); %D表示各个城市之间的距离
g_route=[NC_max,n+20]; %各代最佳路线
temp_best_length=inf.*ones(NC_max,1); %各代最佳路线的长度
length_ave=zeros(NC_max,1); %各代路线的平均长度
temp_pra=0;
%---------------绘制仓库和客户点--------------------%
figure(1)
plot(coordinate(1,1),coordinate(1,2),'r-*');
text(coordinate(1,1),coordinate(1,2),[ ' ''仓库'] )
holdon
fori=2:n
plot(coordinate(i,1),coordinate(i,2),'-*');
text( coordinate(i,1),coordinate(i,2),[' ' num2str(i-1)] )
hold on
end
%--------------------求距离--------------------
fori=1:n
for j=1:n
if i~=j
Distance(i,j)=((coordinate(i,1)-coordinate(j,1))^2+(coordinate(i,2)-coordinate(j,2))^2)^0.5;
else
Distance(i,j)=eps;
end
Distance(j,i)=Distance(i,j);
end
end
% 构造工厂之间比分别从工厂出发时所能节省的路程
U=zeros(n,n);
fori=1:n
for j=1:n
if i~=j
U(i,j)=Distance(i,1)+Distance(j,1)-Distance(i,j);
else
U(i,j)=eps;
end
U(j,i)=U(i,j);
end
end
load=0; %
Eta=1./Distance; %Eta为启发因子，离散模型设为距离的倒数
Tau=ones(n,n); %Tau为信息素矩阵，行：前一个城市的标号；列：后一个城市的标号。
Tabu=zeros(m,n+15); %存储并记录路径的生成
NC=1; %迭代计数器
%------------------迭代开始------------------%
whileNC<=NC_max
Tabu(:,1)=1; %每只蚂蚁都是从仓库（及节点1）开始出发
for i=1:m
visited=Tabu(i,:);
visited=visited(visited>0);
to_visit=setdiff(1:n,visited); %在集合论中，to_visit = (1:n) - visited
c_temp=length(to_visit);
j=1;
%------------------------------------------------------------%
% 第一步按概率选择城市，第二步根据车子的载重量限制判定是否选择这个城市%
%------------------------------------------------------------%
while j<=n %参观了的城市数目
if ~isempty(to_visit) %判断to_visit不是空，是否还有城市没有访问到
%-----------------求解visited的最后一个城市，到所有还没有参观的城市的概率-----------------%
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);
end
temp_pra=rand;
if temp_pra<Prab %如果满足伪随机选择概率，直接选择概率最大的，这样收敛更快
Select=find(max(x));
else
x=x/(sum(x)); %用轮盘赌法选择下一个城市
xcum=cumsum(x);
Select=find(xcum>=rand);
end
if isempty(Select) %如果选不了下一个城市了，就回到1点
Select=1;
load=load+demand(Select); %select=1；t（select）=0
else
load=load+demand(to_visit(Select(1))); %select表示城市标号，to_visit(select)才表示城市
end
end
visited=Tabu(i,:);
visited=visited(visited>0); %已经拜访过的城市
to_visit=setdiff(1:n,visited);
x=[];
end
load=0;
end
%------------------求蚁群m只蚂蚁路径的路程------------------%
L=zeros(m,1); %共有m只蚂蚁
for i=1:m
temp=Tabu(i,:);
R=temp(temp>0);
for j=1:(length(R)-1)
L(i)=L(i)+Distance(R(j),R(j+1));
end
end
minL = min(L);
%%------------------引入模拟退火，产生更新集-------------------%
[Tabu,L]=SA(Tabu,L,m,n,demand,W,Distance,minL,t);
%---------模拟退火引入完成------------%
temp_best_length(NC)=min(L);
pos=find(L==temp_best_length(NC));
g_route(NC,1:length(Tabu(pos(1),:)))=Tabu(pos(1),:);
%% 画图，做成视频
figure(2)
DrawRoute(coordinate,temp);
title(['进化第',num2str(NC),'代']);
PICTURE(NC)=getframe;
hold off
%%2-opt优化完毕
length_ave(NC)=mean(L);
NC=NC+1 %显示进化的代数
%------------------更新信息素-------------------%
Delta_Tau=zeros(n,n);
for i=1:m
temp=Tabu(i,:);
R=temp(temp>0);
for j=1:(length(R)-1)
Delta_Tau(R(j),R(j+1))=Delta_Tau(R(j),R(j+1))+Q/L(i);
end
end
Tau=(1-Rho).*Tau+Delta_Tau;
%----------------更新模拟退后的温度---------------
t=t*down;
if t<tmin
t = tmax;
end
%------------------路径矩阵清零-------------------%
Tabu=zeros(m,n);
load=0;
end
%-------------------输出结果----------------------%
Pos=find(temp_best_length==min(temp_best_length)); %找到最短路径对应的行
best_route=g_route(Pos(1),:); %最短的路径
best_length=temp_best_length(Pos(1)); %最短路径的长度
%画图
best_route=best_route(best_route>0);
disp(['最短路线为：',num2str(best_route)])
disp(['最短路线对应的线路长度为：',num2str(best_length)])
%-------------------绘制迭代优化结果图----------------------%
%figure()
% p= plot(temp_best_length);
%set(p,'Color','red','LineWidth',2);
%xlabel('迭代代数')
%ylabel('目标函数值')
%title('进化过程')
%-------------------绘制迅游过程图----------------------%
figure(3)
plot(coordinate(1,1),coordinate(1,2),'r-*');
text( coordinate(1,1),coordinate(1,2),[' ''仓库'] )
hold on
fori=2:n
plot(coordinate(i,1),coordinate(i,2),'-*');
text( coordinate(i,1),coordinate(i,2),[' ' num2str(i-1)] )
hold on
end
figure(3)
%plot([coordinate(best_route,1)],[coordinate(best_route,2)],'r-*')
DrawRoute(coordinate,best_route);
holdon
PICTURE(NC_max)= getframe;
%制作视频
movie2avi(PICTURE,'SA-ACO-VRP.avi','fps',3,'quality',90)%输出视频
%movie2avi(PICTURE,'SA-ACO.avi', 'compression', 'None')%输出视频
结果：如有疑问，请咨询qq：778961303。
