【VRP问题】基于节约算法求解TWVRP问题

车辆路径问题（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维决策变量

C-W节约算法：

基本思想是把各点单独与货源相连，构成若干条仅含一个配送点的线路，总费用为两倍从原点到各点的距离费用；然后计算将点 i 和点 j 连接在一条线路上费用节约值：

S(i,j)=Coi+Cio+Coj+Cjo−(Coi+Cij+Cjo)=Coi+Coj+CijS(i,j)=Coi+Cio+Coj+Cjo−(Coi+Cij+Cjo)=Coi+Coj+Cij

具体步骤：
（1）计算节约值S(i,j)，按从大到小排序
（2）考虑表格中最大元素 Smax（i,j）Smax（i,j），对应点i和j，按条件进行操作：
1. 若i和j均不在构成线路上，则得到线路 o -> i ->j ->o，转到（3）
2. 若i或j在已构成线路上，但不是内点 0 -> i ->o，则可连接，转到（3）
3. 若i和j位于已构成不同线路上，且均不是内点，则连接得到线路，转到（3）
4. 若i和j位于已构成的同一线路，则不连接，转到（3）
（3）划去第i行和第j列，即i点不能再到其他点，j点也不能由其他店到达
（4）若所有元素均被划去，则得到完整线路，算法终止；否则，在没有划去的元素中选最大元素，转至（2）。

%
 
clear
clc
tic
%% 用importdata这个函数来读取文件 
rc208=importdata('rc208.txt');
cap=1000;
%% 提取数据信息
E=rc208(1,5);                                                    %配送中心时间窗开始时间
L=rc208(1,6);                                                    %配送中心时间窗结束时间
vertexs=rc208(:,2:3);                                            %所有点的坐标x和y
customer=vertexs(2:end,:);                                       %顾客坐标
cusnum=size(customer,1);                                         %顾客数
demands=rc208(2:end,4);                                          %需求量
a=rc208(2:end,5);                                                %顾客时间窗开始时间[a[i],b[i]]
b=rc208(2:end,6);                                                %顾客时间窗结束时间[a[i],b[i]]
s=rc208(2:end,7);                                                %客户点的服务时间
h=pdist(vertexs);
dist=squareform(h);                                             %距离矩阵，满足三角关系，暂用距离表示花费c[i][j]=dist[i][j]
%% CW法构造VRPTW初始解
[init_vc,init_TD,init_vl,violate_INTW]=init_TW(rc208,cap);
initNV=size(init_vc,1);
str1=['车辆行驶总距离 =  ' num2str(init_TD)];
disp(str1)
str2=['车辆使用数目 =  ' num2str(initNV)];
disp(str2)
%% 判断最优解是否满足时间窗约束和载重量约束，0表示违反约束，1表示满足全部约束
flag=Judge(init_vc,cap,demands,a,b,L,s,dist);
%% 检查最优解中是否存在元素丢失的情况，丢失元素，如果没有则为空
DEL=Judge_Del(init_vc);
%% 画出配送路线图
vertexs=rc208(:,2:3);                                            %所有点的坐标x和y
draw_Best(init_vc,vertexs);
toc