6 formulations implemented:
7 - mtz -- Miller-Tucker-Zemlin's potential formulation
8 - mtz_strong -- Miller-Tucker-Zemlin's potential formulation with stronger constraint
9 - scf -- single-commodity flow formulation
10 - mcf -- multi-commodity flow formulation
12 Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
14 from pyscipopt
import Model, quicksum, multidict
17 """mtz: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
18 (potential formulation)
21 - c[i,j]: cost for traversing arc (i,j)
22 Returns a model, ready to be solved.
25 model = Model(
"atsp - mtz")
28 for i
in range(1,n+1):
29 u[i] = model.addVar(lb=0, ub=n-1, vtype=
"C", name=
"u(%s)"%i)
30 for j
in range(1,n+1):
32 x[i,j] = model.addVar(vtype=
"B", name=
"x(%s,%s)"%(i,j))
34 for i
in range(1,n+1):
35 model.addCons(
quicksum(x[i,j]
for j
in range(1,n+1)
if j != i) == 1,
"Out(%s)"%i)
36 model.addCons(
quicksum(x[j,i]
for j
in range(1,n+1)
if j != i) == 1,
"In(%s)"%i)
38 for i
in range(1,n+1):
39 for j
in range(2,n+1):
41 model.addCons(u[i] - u[j] + (n-1)*x[i,j] <= n-2,
"MTZ(%s,%s)"%(i,j))
43 model.setObjective(
quicksum(c[i,j]*x[i,j]
for (i,j)
in x),
"minimize")
51 """mtz_strong: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
52 (potential formulation, adding stronger constraints)
55 c[i,j] - cost for traversing arc (i,j)
56 Returns a model, ready to be solved.
59 model = Model(
"atsp - mtz-strong")
62 for i
in range(1,n+1):
63 u[i] = model.addVar(lb=0, ub=n-1, vtype=
"C", name=
"u(%s)"%i)
64 for j
in range(1,n+1):
66 x[i,j] = model.addVar(vtype=
"B", name=
"x(%s,%s)"%(i,j))
68 for i
in range(1,n+1):
69 model.addCons(
quicksum(x[i,j]
for j
in range(1,n+1)
if j != i) == 1,
"Out(%s)"%i)
70 model.addCons(
quicksum(x[j,i]
for j
in range(1,n+1)
if j != i) == 1,
"In(%s)"%i)
72 for i
in range(1,n+1):
73 for j
in range(2,n+1):
75 model.addCons(u[i] - u[j] + (n-1)*x[i,j] + (n-3)*x[j,i] <= n-2,
"LiftedMTZ(%s,%s)"%(i,j))
77 for i
in range(2,n+1):
78 model.addCons(-x[1,i] - u[i] + (n-3)*x[i,1] <= -2, name=
"LiftedLB(%s)"%i)
79 model.addCons(-x[i,1] + u[i] + (n-3)*x[1,i] <= n-2, name=
"LiftedUB(%s)"%i)
81 model.setObjective(
quicksum(c[i,j]*x[i,j]
for (i,j)
in x),
"minimize")
88 """scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem
91 - c[i,j]: cost for traversing arc (i,j)
92 Returns a model, ready to be solved.
94 model = Model(
"atsp - scf")
97 for i
in range(1,n+1):
98 for j
in range(1,n+1):
100 x[i,j] = model.addVar(vtype=
"B", name=
"x(%s,%s)"%(i,j))
102 f[i,j] = model.addVar(lb=0, ub=n-1, vtype=
"C", name=
"f(%s,%s)"%(i,j))
104 f[i,j] = model.addVar(lb=0, ub=n-2, vtype=
"C", name=
"f(%s,%s)"%(i,j))
106 for i
in range(1,n+1):
107 model.addCons(
quicksum(x[i,j]
for j
in range(1,n+1)
if j != i) == 1,
"Out(%s)"%i)
108 model.addCons(
quicksum(x[j,i]
for j
in range(1,n+1)
if j != i) == 1,
"In(%s)"%i)
110 model.addCons(
quicksum(f[1,j]
for j
in range(2,n+1)) == n-1,
"FlowOut")
112 for i
in range(2,n+1):
113 model.addCons(
quicksum(f[j,i]
for j
in range(1,n+1)
if j != i) - \
114 quicksum(f[i,j]
for j
in range(1,n+1)
if j != i) == 1,
"FlowCons(%s)"%i)
116 for j
in range(2,n+1):
117 model.addCons(f[1,j] <= (n-1)*x[1,j],
"FlowUB(%s,%s)"%(1,j))
118 for i
in range(2,n+1):
120 model.addCons(f[i,j] <= (n-2)*x[i,j],
"FlowUB(%s,%s)"%(i,j))
122 model.setObjective(
quicksum(c[i,j]*x[i,j]
for (i,j)
in x),
"minimize")
130 """mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem
133 - c[i,j]: cost for traversing arc (i,j)
134 Returns a model, ready to be solved.
139 for i
in range(1,n+1):
140 for j
in range(1,n+1):
142 x[i,j] = model.addVar(vtype=
"B", name=
"x(%s,%s)"%(i,j))
143 if i != j
and j != 1:
144 for k
in range(2,n+1):
146 f[i,j,k] = model.addVar(ub=1, vtype=
"C", name=
"f(%s,%s,%s)"%(i,j,k))
148 for i
in range(1,n+1):
149 model.addCons(
quicksum(x[i,j]
for j
in range(1,n+1)
if j != i) == 1,
"Out(%s)"%i)
150 model.addCons(
quicksum(x[j,i]
for j
in range(1,n+1)
if j != i) == 1,
"In(%s)"%i)
152 for k
in range(2,n+1):
153 model.addCons(
quicksum(f[1,i,k]
for i
in range(2,n+1)
if (1,i,k)
in f) == 1,
"FlowOut(%s)"%k)
154 model.addCons(
quicksum(f[i,k,k]
for i
in range(1,n+1)
if (i,k,k)
in f) == 1,
"FlowIn(%s)"%k)
156 for i
in range(2,n+1):
158 model.addCons(
quicksum(f[j,i,k]
for j
in range(1,n+1)
if (j,i,k)
in f) == \
159 quicksum(f[i,j,k]
for j
in range(1,n+1)
if (i,j,k)
in f),
160 "FlowCons(%s,%s)"%(i,k))
163 model.addCons(f[i,j,k] <= x[i,j],
"FlowUB(%s,%s,%s)"%(i,j,k))
165 model.setObjective(
quicksum(c[i,j]*x[i,j]
for (i,j)
in x),
"minimize")
173 """sequence: make a list of cities to visit, from set of arcs"""
179 for i
in range(len(arcs)-2):
185 if __name__ ==
"__main__":
187 c = { (1,1):0, (1,2):1989, (1,3):102, (1,4):102, (1,5):103,
188 (2,1):104, (2,2):0, (2,3):11, (2,4):104, (2,5):108,
189 (3,1):107, (3,2):108, (3,3):0, (3,4):19, (3,5):102,
190 (4,1):109, (4,2):102, (4,3):107, (4,4):0, (4,5):15,
191 (5,1):13, (5,2):103, (5,3):104, (5,4):101, (5,5):0,
197 cost = model.getObjVal()
199 print(
"Miller-Tucker-Zemlin's model:")
200 print(
"Optimal value:", cost)
202 for v
in model.getVars():
203 if model.getVal(v) > 0.001:
204 print(v.name,
"=", model.getVal(v))
207 sol = [i
for (p,i)
in sorted([(int(model.getVal(u[i])+.5),i)
for i
in range(1,n+1)])]
209 arcs = [(i,j)
for (i,j)
in x
if model.getVal(x[i,j]) > .5]
217 cost = model.getObjVal()
219 print(
"Miller-Tucker-Zemlin's model with stronger constraints:")
220 print(
"Optimal value:",cost)
222 for v
in model.getVars():
223 if model.getVal(v) > 0.001:
224 print(v.name,
"=", model.getVal(v))
227 sol = [i
for (p,i)
in sorted([(int(model.getVal(u[i])+.5),i)
for i
in range(1,n+1)])]
229 arcs = [(i,j)
for (i,j)
in x
if model.getVal(x[i,j]) > .5]
237 cost = model.getObjVal()
239 print(
"Single-commodity flow formulation:")
240 print(
"Optimal value:",cost)
242 for v
in model.getVars():
243 if model.getVal(v) > 0.001:
244 print(v.name,
"=", model.getVal(v))
247 arcs = [(i,j)
for (i,j)
in x
if model.getVal(x[i,j]) > .5]
255 cost = model.getObjVal()
257 print(
"Multi-commodity flow formulation:")
258 print(
"Optimal value:",cost)
260 for v
in model.getVars():
261 if model.getVal(v)>0.001:
262 print(v.name,
"=", model.getVal(v))
265 arcs = [(i,j)
for (i,j)
in x
if model.getVal(x[i,j]) > .5]