一、简介
模拟退火算法
著名的模拟退火算法，它是一种基于蒙特卡洛思想设计的近似求解最优化问题的方法。
一点历史——如果你不感兴趣，可以跳过
美国物理学家 N.Metropolis 和同仁在1953年发表研究复杂系统、计算其中能量分布的文章，他们使用蒙特卡罗模拟法计算多分子系统中分子的能量分布。这相当于是本文所探讨之问题的开始，事实上，模拟退火中常常被提到的一个名词就是Metropolis准则，后面我们还会介绍。
美国IBM公司物理学家 S.Kirkpatrick、C. D. Gelatt 和 M. P. Vecchi 于1983年在《Science》上发表了一篇颇具影响力的文章：《以模拟退火法进行最优化（Optimization by Simulated Annealing）》。他们借用了Metropolis等人的方法探讨一种旋转玻璃态系统（spin glass system）时，发觉其物理系统的能量和一些组合最优（combinatorial optimization）问题（著名的旅行推销员问题TSP即是一个代表例子）的成本函数相当类似：寻求最低成本即似寻求最低能量。由此，他们发展出以 Metropolis 方法为本的一套算法，并用其来解决组合问题等的寻求最优解。
几乎同时，欧洲物理学家 V.Carny 也发表了几乎相同的成果，但两者是各自独立发现的；只是Carny“运气不佳”，当时没什么人注意到他的大作；或许可以说，《Science》杂志行销全球，“曝光度”很高，素负盛名，而Carny却在另外一本发行量很小的专门学术期刊《J.Opt.Theory Appl.》发表其成果因而并未引起应有的关注。
Kirkpatrick等人受到Metropolis等人用蒙特卡罗模拟的启发而发明了“模拟退火”这个名词，因为它和物体退火过程相类似。寻找问题的最优解（最值）即类似寻找系统的最低能量。因此系统降温时，能量也逐渐下降，而同样意义地，问题的解也“下降”到最值。
一、什么是退火——物理上的由来
在热力学上，退火（annealing）现象指物体逐渐降温的物理现象，温度愈低，物体的能量状态会低；够低后，液体开始冷凝与结晶，在结晶状态时，系统的能量状态最低。大自然在缓慢降温（亦即，退火）时，可“找到”最低能量状态：结晶。但是，如果过程过急过快，快速降温（亦称「淬炼」，quenching）时，会导致不是最低能态的非晶形。
如下图所示，首先（左图）物体处于非晶体状态。我们将固体加温至充分高（中图），再让其徐徐冷却，也就退火（右图）。加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小（此时物体以晶体形态呈现）。
似乎，大自然知道慢工出细活：缓缓降温，使得物体分子在每一温度时，能够有足够时间找到安顿位置，则逐渐地，到最后可得到最低能态，系统最安稳。
二、模拟退火（Simulate Anneal）
如果你对退火的物理意义还是晕晕的，没关系我们还有更为简单的理解方式。想象一下如果我们现在有下面这样一个函数，现在想求函数的（全局）最优解。如果采用Greedy策略，那么从A点开始试探，如果函数值继续减少，那么试探过程就会继续。而当到达点B时，显然我们的探求过程就结束了（因为无论朝哪个方向努力，结果只会越来越大）。最终我们只能找打一个局部最后解B。
模拟退火其实也是一种Greedy算法，但是它的搜索过程引入了随机因素。模拟退火算法以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。以上图为例，模拟退火算法在搜索到局部最优解B后，会以一定的概率接受向右继续移动。也许经过几次这样的不是局部最优的移动后会到达B 和C之间的峰点，于是就跳出了局部最小值B。
根据Metropolis准则，粒子在温度T时趋于平衡的概率为exp(-ΔE/(kT))，其中E为温度T时的内能，ΔE为其改变数,k为Boltzmann常数。Metropolis准则常表示为
Metropolis准则表明，在温度为T时，出现能量差为dE的降温的概率为P(dE)，表示为：P(dE) = exp( dE/(kT) )。其中k是一个常数，exp表示自然指数，且dE<0。所以P和T正相关。这条公式就表示：温度越高，出现一次能量差为dE的降温的概率就越大；温度越低，则出现降温的概率就越小。又由于dE总是小于0（因为退火的过程是温度逐渐下降的过程），因此dE/kT < 0 ，所以P(dE)的函数取值范围是(0,1) 。随着温度T的降低，P(dE)会逐渐降低。
我们将一次向较差解的移动看做一次温度跳变过程，我们以概率P(dE)来接受这样的移动。也就是说，在用固体退火模拟组合优化问题，将内能E模拟为目标函数值 f，温度T演化成控制参数 t，即得到解组合优化问题的模拟退火演算法：由初始解 i 和控制参数初值 t 开始，对当前解重复“产生新解→计算目标函数差→接受或丢弃”的迭代，并逐步衰减 t 值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling Schedule)控制，包括控制参数的初值 t 及其衰减因子Δt 、每个 t 值时的迭代次数L和停止条件S。
总结起来就是：
若f( Y(i+1) ) <= f( Y(i) ) (即移动后得到更优解)，则总是接受该移动；
若f( Y(i+1) ) > f( Y(i) ) (即移动后的解比当前解要差)，则以一定的概率接受移动，而且这个概率随着时间推移逐渐降低（逐渐降低才能趋向稳定）相当于上图中，从B移向BC之间的小波峰时，每次右移（即接受一个更糟糕值）的概率在逐渐降低。如果这个坡特别长，那么很有可能最终我们并不会翻过这个坡。如果它不太长，这很有可能会翻过它，这取决于衰减 t 值的设定。
关于普通Greedy算法与模拟退火，有一个有趣的比喻：
普通Greedy算法：兔子朝着比现在低的地方跳去。它找到了不远处的最低的山谷。但是这座山谷不一定最低的。这就是普通Greedy算法，它不能保证局部最优值就是全局最优值。
模拟退火：兔子喝醉了。它随机地跳了很长时间。这期间，它可能走向低处，也可能踏入平地。但是，它渐渐清醒了并朝最低的方向跳去。这就是模拟退火。
车辆路径问题（VRP）
车辆路径问题（VRP）最早是由 Dantzig 和 Ramser 于1959年首次提出，它是指一定数量的客户，各自有不同数量的货物需求，配送中心向客户提供货物，由一个车队负责分送货物，组织适当的行车路线，目标是使得客户的需求得到满足，并能在一定的约束下，达到诸如路程最短、成本最小、耗费时间最少等目的。
带时间窗的车辆路径问题(VRPTW)
由于VRP问题的持续发展，考虑需求点对于车辆到达的时间有所要求之下，在车辆途程问题之中加入时窗的限制，便成为带时间窗车辆路径问题（VRP with Time Windows, VRPTW）。带时间窗车辆路径问题（VRPTW）是在VRP上加上了客户的被访问的时间窗约束。在VRPTW问题中，除了行驶成本之外, 成本函数还要包括由于早到某个客户而引起的等待时间和客户需要的服务时间。在VRPTW中，车辆除了要满足VRP问题的限制之外，还必须要满足需求点的时窗限制，而需求点的时窗限制可以分为两种，一种是硬时窗（Hard Time Window），硬时窗要求车辆必须要在时窗内到达，早到必须等待，而迟到则拒收；另一种是软时窗（Soft Time Window），不一定要在时窗内到达，但是在时窗之外到达必须要处罚，以处罚替代等待与拒收是软时窗与硬时窗最大的不同。
模型1问题定义：
The VRPTW is defined by a fleet of vehicles K={1,…,k} , a set of customers C={1,…,n} , and a directed graph G, Typically the fleet is considered to be homogeneous, that is, all vehicles are identical. The graph consists of |C| + 2 vertices, where the customers are denoted 1, 2,…,n and the depot is represented by the vertices 0 (“the starting depot”) and n + 1 (“the returning depot”). The set of all vertices, that is, 0,1,… , n+1 is denoted N. The set of arcs, A, represents direct connections between the depot and the customers and among the customers. There are no arcs ending at vertex 0 or originating from vertex n + 1. With each arc (i,j), where i≠j , we associate a cost cij and a time tij , which may include service time at customer i.
Each vehicle has a capacity Q and each customer i a demand qi . Each customer i has a time window [ai,bi] and a vehicle must arrive at the customer before bi . If it arrives before the time window opens, it has to wait until ai to service the customer. The time windows for both depots are assumed to be identical to [a0,b0] which represents the scheduling horizon. The vehicles may not leave the depot before a0 and must return at the latest at time bn+1 .
It is assumed that Q, ai , bi , qi , cij are non-negative integers and tij are positive integers. Note that this assumption is necessary to develop an algorithm for the shortest path with resource constraints used in the column generation approach presented later. Furthermore it is assumed that the triangle inequality is satisfied for both cij and tij .
The decision variable s ik is defined for each vertex i and each vehicle k and denotes the time vehicle k starts to service customer i. In case vehicle k does not service customer i, s ik has no meaning and consequently it’s value is considered irrelevant. We assume a0 = 0 and therefore s 0k = 0, for all k.
模型2(参考2017 A generalized formulation for vehicle routing problems)：
该模型为2维决策变量
三、源代码
tic
clear
clc
%% 用importdata这个函数来读取文件
c101=importdata('c101.txt');
cap=200;
%% 提取数据信息
E=c101(1,5); %配送中心时间窗开始时间
L=c101(1,6); %配送中心时间窗结束时间
vertexs=c101(:,2:3); %所有点的坐标x和y
customer=vertexs(2:end,:); %顾客坐标
cusnum=size(customer,1); %顾客数
v_num=5; %车辆最多使用数目
demands=c101(2:end,4); %需求量
a=c101(2:end,5); %顾客时间窗开始时间[a[i],b[i]]
b=c101(2:end,6); %顾客时间窗结束时间[a[i],b[i]]
s=c101(2:end,7); %客户点的服务时间
h=pdist(vertexs);
dist=squareform(h); %距离矩阵
%% 模拟退火参数
belta=10; %违反的容量约束的惩罚函数系数
gama=100; %违反时间窗约束的惩罚函数系数
MaxOutIter=2000; %外层循环最大迭代次数
MaxInIter=300; %里层循环最大迭代次数
T0=100; %初始温度
alpha=0.99; %冷却因子
pSwap=0.2; %选择交换结构的概率
pReversion=0.5; %选择逆转结构的概率
pInsertion=1-pSwap-pReversion; %选择插入结构的概率
N=cusnum+v_num-1; %解长度=顾客数目+车辆最多使用数目-1
%% 随机构造初始解
currS=randperm(N); %随机构造初始解
[currVC,NV,TD,violate_num,violate_cus]=decode(currS,cusnum,cap,demands,a,b,L,s,dist); %对初始解解码
%求初始配送方案的成本=车辆行驶总成本+belta*违反的容量约束之和+gama*违反时间窗约束之和
currCost=costFuction(currVC,a,b,s,L,dist,demands,cap,belta,gama);
Sbest=currS; %初始将全局最优解赋值为初始解
bestVC=currVC; %初始将全局最优配送方案赋值为初始配送方案
bestCost=currCost; %初始将全局最优解的总成本赋值为初始解总成本
BestCost=zeros(MaxOutIter,1); %记录每一代全局最优解的总成本
T=T0; %温度初始化
%% 模拟退火
for outIter=1:MaxOutIter
for inIter=1:MaxInIter
newS=Neighbor(currS,pSwap,pReversion,pInsertion); %经过邻域结构后产生的新的解
%对新解解码=车辆行驶总成本+belta*违反的容量约束之和+gama*违反时间窗约束之和
newVC=decode(newS,cusnum,cap,demands,a,b,L,s,dist);
newCost=costFuction(newVC,a,b,s,L,dist,demands,cap,belta,gama); %求初始配送方案的成本
%如果新解比当前解更好，则更新当前解，以及当前解的总成本
if newCost<=currCost
currS=newS;
currVC=newVC;
currCost=newCost;
else
%如果新解不如当前解好，则采用退火准则，以一定概率接受新解
delta=(newCost-currCost)/currCost; %计算新解与当前解总成本相差的百分比
P=exp(-delta/T); %计算接受新解的概率
%如果0~1的随机数小于P，则接受新解，并更新当前解，以及当前解总成本
if rand<=P
currS=newS;
currVC=newVC;
currCost=newCost;
end
end
%将当前解与全局最优解进行比较，如果当前解更好，则更新全局最优解，以及全局最优解总成本
if currCost<=bestCost
Sbest=currS;
bestVC=currVC;
bestCost=currCost;
end
end
%记录外层循环每次迭代的全局最优解的总成本
BestCost(outIter)=bestCost;
%显示外层循环每次迭代的信全局最优解的总成本
disp(['第',num2str(outIter),'代全局最优解:'])
[bestVC,bestNV,bestTD,best_vionum,best_viocus]=decode(Sbest,cusnum,cap,demands,a,b,L,s,dist); %对全局最优解解码
disp(['车辆使用数目：',num2str(bestNV),'，车辆行驶总距离：',num2str(bestTD),'，违反约束路径数目：',num2str(best_vionum),'，违反约束顾客数目：',num2str(best_viocus)]);
fprintf('\n')
%更新当前温度
T=alpha*T;
end
%% 打印外层循环每次迭代的全局最优解的总成本变化趋势图
figure;
plot(BestCost,'LineWidth',1);
title('全局最优解的总成本变化趋势图')
xlabel('迭代次数');
ylabel('总成本');
%% 打印全局最优解路线图
draw_Best(bestVC,vertexs);
toc
三、运行结果
四、备注
版本：2014a
