-
Notifications
You must be signed in to change notification settings - Fork 253
/
vrp.py
111 lines (94 loc) · 3.38 KB
/
vrp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
"""
vrp.py: solve the vehicle routing problem.
approach:
- start with assignment model
- add cuts until all components of the graph are connected
Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
"""
import math
import random
import networkx
from pyscipopt import Model, quicksum, multidict
def solve_vrp(V,c,m,q,Q):
"""solve_vrp -- solve the vehicle routing problem.
- start with assignment model (depot has a special status)
- add cuts until all components of the graph are connected
Parameters:
- V: set/list of nodes in the graph
- c[i,j]: cost for traversing edge (i,j)
- m: number of vehicles available
- q[i]: demand for customer i
- Q: vehicle capacity
Returns the optimum objective value and the list of edges used.
"""
def addcut(cut_edges):
"""addcut: add constraint to eliminate infeasible solutions
Parameters:
- cut_edges: list of edges in the current solution, except connections to depot
Returns True if a cut was added, False otherwise
"""
G = networkx.Graph()
G.add_edges_from(cut_edges)
Components = networkx.connected_components(G)
cut = False
for S in Components:
S_card = len(S)
q_sum = sum(q[i] for i in S)
NS = int(math.ceil(float(q_sum)/Q))
S_edges = [(i,j) for i in S for j in S if i<j and (i,j) in cut_edges]
if S_card >= 3 and (len(S_edges) >= S_card or NS > 1):
add = model.addCons(quicksum(x[i,j] for i in S for j in S if j > i) <= S_card-NS)
cut = True
return cut
model = Model("vrp")
x = {}
for i in V:
for j in V:
if j > i and i == V[0]: # depot
x[i,j] = model.addVar(ub=2, vtype="I", name="x(%s,%s)"%(i,j))
elif j > i:
x[i,j] = model.addVar(ub=1, vtype="I", name="x(%s,%s)"%(i,j))
model.addCons(quicksum(x[V[0],j] for j in V[1:]) == 2*m, "DegreeDepot")
for i in V[1:]:
model.addCons(quicksum(x[j,i] for j in V if j < i) +
quicksum(x[i,j] for j in V if j > i) == 2, "Degree(%s)"%i)
model.setObjective(quicksum(c[i,j]*x[i,j] for i in V for j in V if j>i), "minimize")
model.hideOutput()
EPS = 1.e-6
while True:
model.optimize()
edges = []
for (i,j) in x:
if model.getVal(x[i,j]) > EPS:
if i != V[0] and j != V[0]:
edges.append((i,j))
if addcut(edges) == False:
break
return model.getObjVal(),edges
def distance(x1,y1,x2,y2):
"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
return math.sqrt((x2-x1)**2 + (y2-y1)**2)
def make_data(n):
"""make_data: compute matrix distance based on euclidean distance"""
V = range(1,n+1)
x = dict([(i,random.random()) for i in V])
y = dict([(i,random.random()) for i in V])
c,q = {},{}
Q = 100
for i in V:
q[i] = random.randint(10,20)
for j in V:
if j > i:
c[i,j] = distance(x[i],y[i],x[j],y[j])
return V,c,q,Q
if __name__ == "__main__":
import sys
n = 19
m = 3
seed = 1
random.seed(seed)
V,c,q,Q = make_data(n)
z,edges = solve_vrp(V,c,m,q,Q)
print("Optimal solution:",z)
print("Edges in the solution:")
print(sorted(edges))