from __future__ import print_function from __future__ import division, print_function from gurobipy import * import re; import math; import matplotlib.pyplot as plt import numpy as np import pandas as pd import copy from matplotlib.lines import lineStyles import time
starttime = time.time()
# function to read data from .txt files defreadData(path, nodeNum): nodeNum = nodeNum; cor_X = [] cor_Y = []
f = open(path, 'r'); lines = f.readlines(); count = 0; # read the info for line in lines: count = count + 1; if(count >= 10and count <= 10 + nodeNum): line = line[:-1] str = re.split(r" +", line) cor_X.append(float(str[2])) cor_Y.append(float(str[3]))
# compute the distance matrix disMatrix = [([0] * nodeNum) for p in range(nodeNum)]; # 初始化距离矩阵的维度,防止浅拷贝 # data.disMatrix = [[0] * nodeNum] * nodeNum]; 这个是浅拷贝,容易重复 for i in range(0, nodeNum): for j in range(0, nodeNum): temp = (cor_X[i] - cor_X[j])**2 + (cor_Y[i] - cor_Y[j])**2; disMatrix[i][j] = (int)(math.sqrt(temp)); # disMatrix[i][j] = 0.1 * (int)(10 * math.sqrt(temp)); # if(i == j): # data.disMatrix[i][j] = 0; # print("%6.0f" % (math.sqrt(temp)), end = " "); temp = 0;
return disMatrix;
defprintData(disMatrix): print("-------cost matrix-------\n"); for i in range(len(disMatrix)): for j in range(len(disMatrix)): #print("%d %d" % (i, j)); print("%6.1f" % (disMatrix[i][j]), end = " "); # print(disMatrix[i][j], end = " "); print();
defreportMIP(model, Routes): if model.status == GRB.OPTIMAL: print("Best MIP Solution: ", model.objVal, "\n") var = model.getVars() for i in range(model.numVars): if(var[i].x > 0): print(var[i].varName, " = ", var[i].x) print("Optimal route:", Routes[i])
defgetValue(var_dict, nodeNum): x_value = np.zeros([nodeNum + 1, nodeNum + 1]) for key in var_dict.keys(): a = key[0] b = key[1] x_value[a][b] = var_dict[key].x
# Callback - use lazy constraints to eliminate sub-tours
# Callback - use lazy constraints to eliminate sub-tours
defsubtourelim(model, where): if(where == GRB.Callback.MIPSOL): # make a list of edges selected in the solution print('model._vars', model._vars) # vals = model.cbGetSolution(model._vars) x_value = np.zeros([nodeNum + 1, nodeNum + 1]) for m in model.getVars(): if(m.varName.startswith('x')): # print(var[i].varName) # print(var[i].varName.split('_')) a = (int)(m.varName.split('_')[1]) b = (int)(m.varName.split('_')[2]) x_value[a][b] = model.cbGetSolution(m) print("solution = ", x_value) # print('key = ', model._vars.keys()) # selected = [] # for i in range(nodeNum): # for j in range(nodeNum): # if(i != j and x_value[i][j] > 0.5): # selected.append((i, j)) # selected = tuplelist(selected) # # selected = tuplelist((i,j) for i in range(nodeNum), for if x_value[i][j] > 0.5) # print('selected = ', selected) # find the shortest cycle in the selected edge list tour = subtour(x_value) print('tour = ', tour) if(len(tour) < nodeNum + 1): # add subtour elimination constraint for every pair of cities in tour print("---add sub tour elimination constraint--") # model.cbLazy(quicksum(model._vars[i][j] # for i in tour # for j in tour # if i != j) # <= len(tour)-1) # LinExpr = quicksum(model._vars[i][j] # for i in tour # for j in tour # if i != j) for i,j in itertools.combinations(tour, 2): print(i,j)
model.cbLazy(quicksum(model._vars[i, j] for i,j in itertools.combinations(tour, 2)) <= len(tour)-1) LinExpr = quicksum(model._vars[i, j] for i,j in itertools.combinations(tour, 2)) print('LinExpr = ', LinExpr) print('RHS = ', len(tour)-1)
# compute the degree of each node in given graph defcomputeDegree(graph): degree = np.zeros(len(graph)) for i in range(len(graph)): for j in range(len(graph)): if(graph[i][j] > 0.5): degree[i] = degree[i] + 1 degree[j] = degree[j] + 1 print('degree', degree) return degree
# given a graph, get the edges of this graph deffindEdges(graph): edges = [] for i in range(1, len(graph)): for j in range(1, len(graph)): if(graph[i][j] > 0.5): edges.append((i, j))
return edges
# Given a tuplelist of edges, find the shortest subtour defsubtour(graph): # compute degree of each node degree = computeDegree(graph) unvisited = [] for i in range(1, len(degree)): if(degree[i] >= 2): unvisited.append(i) cycle = range(0, nodeNum + 1) # initial length has 1 more city
edges = findEdges(graph) edges = tuplelist(edges) print(edges) while unvisited: # true if list is non-empty thiscycle = [] neighbors = unvisited while neighbors: # true if neighbors is non-empty current = neighbors[0] thiscycle.append(current) unvisited.remove(current) neighbors = [j for i,j in edges.select(current,'*') if j in unvisited] neighbors2 = [i for i,j in edges.select('*',current) if i in unvisited] if(neighbors2): neighbors.extend(neighbors2) # print('current:', current, '\n neighbors', neighbors)
isLink = ((thiscycle[0], thiscycle[-1]) in edges) or ((thiscycle[-1], thiscycle[0]) in edges) if(len(cycle) > len(thiscycle) and len(thiscycle) >= 3and isLink): # print('in = ', ((thiscycle[0], thiscycle[-1]) in edges) or ((thiscycle[-1], thiscycle[0]) in edges)) cycle = thiscycle return cycle return cycle
然后是建模部分的代码,建模部分相比学运筹的人比较熟悉,这里比较特殊的就是求解时候的几行代码:
model.Params.lazyConstraints = 1 : set lazy constraints Parameter
model.optimize(subtourelim) : use callback function when executing branch and bound algorithm
from __future__ import print_function from __future__ import division, print_function from gurobipy import * import re; import math; import matplotlib.pyplot as plt import numpy as np import pandas as pd import copy from matplotlib.lines import lineStyles import time
starttime = time.time()
# function to read data from .txt files defreadData(path, nodeNum): nodeNum = nodeNum; cor_X = [] cor_Y = []
f = open(path, 'r'); lines = f.readlines(); count = 0; # read the info for line in lines: count = count + 1; if(count >= 10and count <= 10 + nodeNum): line = line[:-1] str = re.split(r" +", line) cor_X.append(float(str[2])) cor_Y.append(float(str[3]))
# compute the distance matrix disMatrix = [([0] * nodeNum) for p in range(nodeNum)]; # 初始化距离矩阵的维度,防止浅拷贝 # data.disMatrix = [[0] * nodeNum] * nodeNum]; 这个是浅拷贝,容易重复 for i in range(0, nodeNum): for j in range(0, nodeNum): temp = (cor_X[i] - cor_X[j])**2 + (cor_Y[i] - cor_Y[j])**2; disMatrix[i][j] = (int)(math.sqrt(temp)); # disMatrix[i][j] = 0.1 * (int)(10 * math.sqrt(temp)); # if(i == j): # data.disMatrix[i][j] = 0; # print("%6.0f" % (math.sqrt(temp)), end = " "); temp = 0;
return disMatrix;
defprintData(disMatrix): print("-------cost matrix-------\n"); for i in range(len(disMatrix)): for j in range(len(disMatrix)): #print("%d %d" % (i, j)); print("%6.1f" % (disMatrix[i][j]), end = " "); # print(disMatrix[i][j], end = " "); print();
defreportMIP(model, Routes): if model.status == GRB.OPTIMAL: print("Best MIP Solution: ", model.objVal, "\n") var = model.getVars() for i in range(model.numVars): if(var[i].x > 0): print(var[i].varName, " = ", var[i].x) print("Optimal route:", Routes[i])
defgetValue(var_dict, nodeNum): x_value = np.zeros([nodeNum + 1, nodeNum + 1]) for key in var_dict.keys(): a = key[0] b = key[1] x_value[a][b] = var_dict[key].x
国内运筹学科普好文还是不太多见,很多都是从1到100的文章。分析讲述一些论文的idea什么的,都是为基础非常好的优秀者们看的。详细的讲述从0到1,如何把基础的东西吃透的文章比较少,让我们这些还不够强的孩子着实举步维艰,听讲座听得懂idea,但是做起来却啥啥也不行。希望以后有干货的,实用的文章越来越多,帮助众多学子解决基本的,底层的疑惑文章。哎,感觉国内OR界博士生们对branch and cut, branch and price, branch and cut and price, benders, DW decomposition等都理解比较深,貌似国内OR人均水平为随手实现上述一系列精确算法如探囊取物。我还是继续好好迎头赶上,抓紧修炼自己吧。希望以后能够真正能够慢慢提高理论和实战水准,荷枪实弹学到真技术,为做有质量的学术打好基础。同时,也希望国内运筹发展越来越好吧。
参考文献
[1]: Desrochers, M., Desrosiers, J., & Solomon, M. (1992). A new optimization algorithm for the vehicle routing problem with time windows. Operations research, 40(2), 342-354. https://doi.org/10.1287/opre.40.2.342