数学规划求解器
性能测试之VRPTW
相比于各种各样的算法，用数学规划求解器求解一些模型可以说是非常简单而有效了。随着CLPEX、Gurobi等各种求解器的出现和求解性能的不断提升，它们在一定程度上已经成为了部分企业乃至学者的偏爱。
但是，求解器真的有这么厉害吗？
小编认为，求解器还是存在着明显的局限性的。
01
问题
要了解VRPTW，我们先来聊聊它的前身——VRP问题
1、什么是VRP
车辆路径问题（Vehicle Routing Problem，VRP）最早是由 Dantzig 和 Ramser 于1959年首次提出，它是指一定数量的客户，各自有不同数量的货物需求，配送中心向客户提供货物，由一个车队负责分送货物，组织适当的行车路线，目标是使得客户的需求得到满足，并能在一定的约束下，达到诸如路程最短、成本最小、耗费时间最少等目的。
由此定义不能看出，旅行商问题（Traveling Saleman Problem，TSP）是VRP的特例，由于Gaery已证明TSP是NP难题，因此VRP也属于NP难题。
在基本车辆路线问题(VRP)的基础上，车辆路线问题在学术研究和实际应用上产生了许多不同的延伸和变化型态，包括时窗限制车辆路线问题(vehicle routing problems with time windows, VRPTW)、追求最佳服务时间的车辆路线问题(VRPDT)、多车种车辆路线问运题(fleet size and mix vehicle routing problems, FSVRP)、车辆多次使用的车辆路线问题(vehicle routing problems with multiple use of vehicle， VRPM)、考虑收集的车辆路线问题(（vehice routing problems with backhauls, VRPB)、随机需求车辆路线问题(vehicle routng problem with stochastic demand,VRPSD)等。
2、什么是VRPTW？
VRPTW，顾名思义，即在VRP的基础上加上了时间窗。
由于VRP问题的持续发展，考虑需求点对于车辆到达的时间有所要求之下，在车辆途程问题之中加入时窗的限制，便成为带时间窗车辆路径问题（VRP with Time Windows, VRPTW）。带时间窗车辆路径问题（VRPTW）是在VRP上加上了客户的被访问的时间窗约束。在VRPTW问题中，除了行驶成本之外, 成本函数还要包括由于早到某个客户而引起的等待时间和客户需要的服务时间。
在VRPTW中，车辆除了要满足VRP问题的限制之外，还必须要满足需求点的时窗限制，而需求点的时窗限制可以分为两种，一种是硬时窗（Hard Time Window），硬时窗要求车辆必须要在时窗内到达，早到必须等待，而迟到则拒收；另一种是软时窗（Soft Time Window），不一定要在时窗内到达，但是在时窗之外到达必须要处罚，以处罚替代等待与拒收是软时窗与硬时窗最大的不同。
模型
02
VRPTW的模型在我们公众号的之前的推文中已有详细介绍，因此不再赘述。在此贴上一篇VRPTW的相关推文，感兴趣的小伙伴可以自行浏览~
干货|十分钟快速掌握CPLEX求解VRPTW数学模型（附JAVA代码及CPLEX安装流程）
03
代码
``定义需要用到的参数：
代码语言：javascript复制class Data:
customerNum = 0
vehicleNum = 0
capacity = 0
cor_X = []
cor_Y = []
demand = []
serviceTime = []
readyTime = []
dueTime = []
disMatrix = [[]] `读取数据，计算距离
代码语言：javascript复制def readData(data, path, customerNum):
data.customerNum = customerNum
f = open(path, 'r')
lines = f.readlines()
count = 0
line = lines[4];
line = line[:-1].strip()
str = re.split(r"s{1,}", line)
data.vehicleNum = int(str[0])
data.capacity = float(str[1])
for i in range(9,len(lines)):
line = lines[i];
count = count 1
line = line[:-1].strip()
str = re.split(r"s{1,}", line)
data.cor_X.append(float(str[1]))
data.cor_Y.append(float(str[2]))
data.demand.append(float(str[3]))
data.readyTime.append(float(str[4]))
data.dueTime.append(float(str[5]))
data.serviceTime.append(float(str[6]))
data.cor_X.append(data.cor_X[0])
data.cor_Y.append(data.cor_Y[0])
data.demand.append(data.demand[0])
data.readyTime.append(data.readyTime[0])
data.dueTime.append(data.dueTime[0])
data.serviceTime.append(data.serviceTime[0])
data.disMatrix = [([0] * (data.customerNum 2)) for p in range(data.customerNum 2)]
for i in range(0, data.customerNum 2):
for j in range(0, data.customerNum 2):
if(i != j):
temp = (data.cor_X[i] - data.cor_X[j]) ** 2 (data.cor_Y[i] - data.cor_Y[j]) ** 2
data.disMatrix[i][j] = math.sqrt(temp)
temp = 0
else:
data.disMatrix[i][j] = math.sqrt(100000.0)
return data读取数据，并定义一些参数：
代码语言：javascript复制data = Data()
path = 'c101.txt' #读取算例数据集
customerNum = 100 #设置客户数量
readData(data, path, customerNum)
BigM = 100000 调用gurobi进行模型的建立与求解：
代码语言：javascript复制x = {} #存放决策变量x_ijk
s = {} #s_ik表示车辆k开始服务客户i的时间
model = Model()
#定义决策变量，并加入模型当中：
for i in range(0,data.customerNum 2):
for k in range(0,data.vehicleNum):
name = 's_' str(i) '_' str(k)
s[i,k] = model.addVar(0, 1500, vtype= GRB.CONTINUOUS, name= name)
for j in range(0,data.customerNum 2):
if(i != j):
name = 'x_' str(i) '_' str(j) '_' str(k)
x[i,j,k] = model.addVar(0, 1, vtype= GRB.BINARY, name= name)
#根据式(1)定义目标函数，并加入模型当中：
obj = LinExpr(0)
for i in range(data.customerNum 2):
for k in range(data.vehicleNum):
for j in range(data.customerNum 2):
if(i != j):
obj.addTerms(data.disMatrix[i][j], x[i,j,k])
model.setObjective(obj, GRB.MINIMIZE)
#约束二：
for i in range(1, data.customerNum 2 - 1):
lhs = LinExpr(0)
for k in range(data.vehicleNum):
for j in range(1, data.customerNum 2):
if(i != j):
lhs.addTerms(1, x[i,j,k])
model.addConstr(lhs == 1, name= 'visit_' str(i))
#约束三：
for k in range(0,data.vehicleNum):
lhs = LinExpr(0)
for j in range(1,data.customerNum 2):
lhs.addTerms(1, x[0,j,k])
model.addConstr(lhs == 1, name= 'vehicle_' str(k))
#约束四：
for k in range(data.vehicleNum):
for h in range(1, data.customerNum 2 - 1):
expr1 = LinExpr(0)
expr2 = LinExpr(0)
for i in range(data.customerNum 2 - 1):
if (h != i):
expr1.addTerms(1, x[i,h,k])
for j in range(1,data.customerNum 2):
if (h != j):
expr2.addTerms(1, x[h,j,k])
model.addConstr(expr1 == expr2, name= 'flow_conservation_' str(i))
expr1.clear()
expr2.clear()
#约束五：
for k in range(data.vehicleNum):
lhs = LinExpr(0)
for j in range(data.customerNum 2 - 1):
lhs.addTerms(1, x[j, data.customerNum 2 - 1, k])
model.addConstr(lhs == 1, name= 'enter_' str(k))
#约束六
for k in range(data.vehicleNum):
for i in range(data.customerNum 2):
for j in range(data.customerNum 2):
if(i != j):
model.addConstr(s[i,k] data.disMatrix[i][j] data.serviceTime[i] - s[j,k]- 2 * BigM 2 * BigM * x[i,j,k] <= 0 , name= 'time_windows_')
#约束七：
for k in range(data.vehicleNum):
for i in range(1,data.customerNum 2 - 1):
lhs1 = LinExpr()
lhs2 = LinExpr()
for j in range(1,data.customerNum 2):
if(i != j):
model.addConstr(data.readyTime[i] * x[i,j,k] <= s[i,k], name= 'ready_time_' str(i))
model.addConstr(data.dueTime[i] (BigM - BigM * x[i,j,k]) >= s[i,k], name= 'due_time_' str(i))
#约束八：
for k in range(data.vehicleNum):
model.addConstr(data.readyTime[0] <= s[0,k], name= 'ready_time_' str(0))
model.addConstr(data.dueTime[0] >= s[0,k], name= 'due_time_' str(0))
model.addConstr(data.readyTime[data.customerNum 1] <= s[data.customerNum 1,k], name= 'ready_time_' str(data.customerNum 1))
model.addConstr(data.dueTime[data.customerNum 1] >= s[data.customerNum 1,k], name= 'due_time_' str(data.customerNum 1))
#约束九
for k in range(data.vehicleNum):
lhs = LinExpr(0)
for i in range(1, data.customerNum 2 - 1):
for j in range(1,data.customerNum 2):
if(i != j):
lhs.addTerms(data.demand[i], x[i,j,k])
model.addConstr(lhs <= data.capacity, name= 'capacity_' str(k))求解模型（solve）并输出解：
代码语言：javascript复制model.optimize()
print("nn-----optimal value-----")
print(model.ObjVal)
print("routes = ")
for k in range(data.vehicleNum):
sys.stdout.write("[")
i = 0;
while i <= customerNum 2 - 1:
if(i == customerNum 1):
sys.stdout.write(str(i))
break
for j in range(1,data.customerNum 2):
if(i != j and x[i,j,k].x > 0):
sys.stdout.write(str(i) ",")
break;
i = j
sys.stdout.write("]")
print("")算例及结果展示
04
所有实验使用Python2.7编程实现，Gurobi版本为9.1，使用服务器信息如下：
算例演示
（需要说明的是，为了使实验的对比效果更加明显，我们以10个点为一组增加数据规模。因此算例二、三、四是由homberger标准算例中的homberger_200_customer_instances截取而来。）
算例一（Solomon标准算例）
输入文件格式为：
结果展示：
算例二（homberger标准算例）
输入文件格式为：
结果展示：
算例三（homberger标准算例）
算例四（homberger标准算例）
输入文件格式为：
结果展示：
（算例下载地址：https://www.sintef.no/projectweb/top/vrptw/homberger-benchmark/1000-customers/）
在运行130个点的算例时，还没有求解完毕代码运行时间就超出了2小时，小编没有再继续运行下去。可以发现，从100个点到130个点仅仅增加了30个点，求解所需时间却在指数级增长。
Gurobi在两个小时内能成功求解的算例规模只有120-130个点，并没有我们想象中的那么大。在企业应用中，更大规模的VRPTW并不少见，但其求解所需时间却不能在企业所能忍受的范围内。
此外，VRPTW其实还算是一个比较简单的路径规划问题，还有很多其他的路径优化问题及其变种，它们比VRPTW更加复杂，如果用Gurobi进行求解，在两个小时内很难达到100个点的数据规模，可能在求解40-50个点的算例时就已经超出时间了。对小规模算例即使求解器可以求解，所需要的求解时间也很长，很多实际场景往往需要几分钟出结果，这个时候求解器就无能为力了。
这足以说明求解器的作用是非常有限的，主要体现在两个方面：
1、其对很多问题无法求解；
2、在理想的时间内能够求解的问题规模不大 。
因此，在解决部分简单且小规模问题时，运用求解器这种相对简单的求解方法当然是上上策，但是在其他大规模问题的求解上，目前的求解器的性能可能还无法在理想时间内实现求解。
END
代码和算例会在留言区给出
