1
1
##@file atsp.py
2
- #@brief solve the asymmetric traveling salesman problem
2
+ # @brief solve the asymmetric traveling salesman problem
3
3
4
4
"""
5
5
11
11
12
12
Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
13
13
"""
14
- from pyscipopt import Model , quicksum , multidict
14
+ from pyscipopt import Model , quicksum
15
15
16
- def mtz (n ,c ):
16
+
17
+ def mtz (n , c ):
17
18
"""mtz: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
18
19
(potential formulation)
19
20
Parameters:
@@ -24,30 +25,29 @@ def mtz(n,c):
24
25
25
26
model = Model ("atsp - mtz" )
26
27
27
- x ,u = {},{}
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 ):
28
+ x , u = {}, {}
29
+ for i in range (1 , n + 1 ):
30
+ u [i ] = model .addVar (lb = 0 , ub = n - 1 , vtype = "C" , name = "u(%s)" % i )
31
+ for j in range (1 , n + 1 ):
31
32
if i != j :
32
- x [i ,j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i ,j ))
33
+ x [i , j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i , j ))
33
34
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 )
35
+ for i in range (1 , n + 1 ):
36
+ model .addCons (quicksum (x [i , j ] for j in range (1 , n + 1 ) if j != i ) == 1 , "Out(%s)" % i )
37
+ model .addCons (quicksum (x [j , i ] for j in range (1 , n + 1 ) if j != i ) == 1 , "In(%s)" % i )
37
38
38
- for i in range (1 ,n + 1 ):
39
- for j in range (2 ,n + 1 ):
39
+ for i in range (1 , n + 1 ):
40
+ for j in range (2 , n + 1 ):
40
41
if i != j :
41
- model .addCons (u [i ] - u [j ] + (n - 1 ) * x [i ,j ] <= n - 2 , "MTZ(%s,%s)" % (i ,j ))
42
+ model .addCons (u [i ] - u [j ] + (n - 1 ) * x [i , j ] <= n - 2 , "MTZ(%s,%s)" % (i , j ))
42
43
43
- model .setObjective (quicksum (c [i ,j ]* x [i ,j ] for (i ,j ) in x ), "minimize" )
44
-
45
- model .data = x ,u
46
- return model
44
+ model .setObjective (quicksum (c [i , j ] * x [i , j ] for (i , j ) in x ), "minimize" )
47
45
46
+ model .data = x , u
47
+ return model
48
48
49
49
50
- def mtz_strong (n ,c ):
50
+ def mtz_strong (n , c ):
51
51
"""mtz_strong: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
52
52
(potential formulation, adding stronger constraints)
53
53
Parameters:
@@ -57,34 +57,34 @@ def mtz_strong(n,c):
57
57
"""
58
58
59
59
model = Model ("atsp - mtz-strong" )
60
-
61
- x ,u = {},{}
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 ):
60
+
61
+ x , u = {}, {}
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 ):
65
65
if i != j :
66
- x [i ,j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i ,j ))
66
+ x [i , j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i , j ))
67
67
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 )
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 )
71
71
72
- for i in range (1 ,n + 1 ):
73
- for j in range (2 ,n + 1 ):
72
+ for i in range (1 , n + 1 ):
73
+ for j in range (2 , n + 1 ):
74
74
if i != j :
75
- model .addCons (u [i ] - u [j ] + (n - 1 )* x [i ,j ] + (n - 3 )* x [j ,i ] <= n - 2 , "LiftedMTZ(%s,%s)" % (i ,j ))
75
+ model .addCons (u [i ] - u [j ] + (n - 1 ) * x [i , j ] + (n - 3 ) * x [j , i ] <= n - 2 , "LiftedMTZ(%s,%s)" % (i , j ))
76
+
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 )
76
80
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" )
80
82
81
- model .setObjective (quicksum (c [i ,j ]* x [i ,j ] for (i ,j ) in x ), "minimize" )
82
-
83
- model .data = x ,u
83
+ model .data = x , u
84
84
return model
85
85
86
86
87
- def scf (n ,c ):
87
+ def scf (n , c ):
88
88
"""scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem
89
89
Parameters:
90
90
- n: number of nodes
@@ -93,40 +93,39 @@ def scf(n,c):
93
93
"""
94
94
model = Model ("atsp - scf" )
95
95
96
- x ,f = {},{}
97
- for i in range (1 ,n + 1 ):
98
- for j in range (1 ,n + 1 ):
96
+ x , f = {}, {}
97
+ for i in range (1 , n + 1 ):
98
+ for j in range (1 , n + 1 ):
99
99
if i != j :
100
- x [i ,j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i ,j ))
101
- if i == 1 :
102
- f [i ,j ] = model .addVar (lb = 0 , ub = n - 1 , vtype = "C" , name = "f(%s,%s)" % (i ,j ))
100
+ x [i , j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i , j ))
101
+ if i == 1 :
102
+ f [i , j ] = model .addVar (lb = 0 , ub = n - 1 , vtype = "C" , name = "f(%s,%s)" % (i , j ))
103
103
else :
104
- f [i ,j ] = model .addVar (lb = 0 , ub = n - 2 , 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 ))
105
105
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 )
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 )
109
109
110
- model .addCons (quicksum (f [1 ,j ] for j in range (2 ,n + 1 )) == n - 1 , "FlowOut" )
110
+ model .addCons (quicksum (f [1 , j ] for j in range (2 , n + 1 )) == n - 1 , "FlowOut" )
111
111
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 )
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 )
115
115
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 ):
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 ):
119
119
if i != j :
120
- model .addCons (f [i ,j ] <= (n - 2 ) * x [i ,j ], "FlowUB(%s,%s)" % (i ,j ))
120
+ model .addCons (f [i , j ] <= (n - 2 ) * x [i , j ], "FlowUB(%s,%s)" % (i , j ))
121
121
122
- model .setObjective (quicksum (c [i ,j ] * x [i ,j ] for (i ,j ) in x ), "minimize" )
122
+ model .setObjective (quicksum (c [i , j ] * x [i , j ] for (i , j ) in x ), "minimize" )
123
123
124
- model .data = x ,f
124
+ model .data = x , f
125
125
return model
126
126
127
127
128
-
129
- def mcf (n ,c ):
128
+ def mcf (n , c ):
130
129
"""mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem
131
130
Parameters:
132
131
- n: number of nodes
@@ -135,134 +134,133 @@ def mcf(n,c):
135
134
"""
136
135
model = Model ("mcf" )
137
136
138
- x ,f = {},{}
139
- for i in range (1 ,n + 1 ):
140
- for j in range (1 ,n + 1 ):
137
+ x , f = {}, {}
138
+ for i in range (1 , n + 1 ):
139
+ for j in range (1 , n + 1 ):
141
140
if i != j :
142
- x [i ,j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i ,j ))
141
+ x [i , j ] = model .addVar (vtype = "B" , name = "x(%s,%s)" % (i , j ))
143
142
if i != j and j != 1 :
144
- for k in range (2 ,n + 1 ):
143
+ for k in range (2 , n + 1 ):
145
144
if i != k :
146
- f [i ,j , k ] = model .addVar (ub = 1 , vtype = "C" , name = "f(%s,%s,%s)" % (i ,j , k ))
145
+ f [i , j , k ] = model .addVar (ub = 1 , vtype = "C" , name = "f(%s,%s,%s)" % (i , j , k ))
147
146
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 )
147
+ for i in range (1 , n + 1 ):
148
+ model .addCons (quicksum (x [i , j ] for j in range (1 , n + 1 ) if j != i ) == 1 , "Out(%s)" % i )
149
+ model .addCons (quicksum (x [j , i ] for j in range (1 , n + 1 ) if j != i ) == 1 , "In(%s)" % i )
151
150
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 )
151
+ for k in range (2 , n + 1 ):
152
+ model .addCons (quicksum (f [1 , i , k ] for i in range (2 , n + 1 ) if (1 , i , k ) in f ) == 1 , "FlowOut(%s)" % k )
153
+ model .addCons (quicksum (f [i , k , k ] for i in range (1 , n + 1 ) if (i , k , k ) in f ) == 1 , "FlowIn(%s)" % k )
155
154
156
- for i in range (2 ,n + 1 ):
155
+ for i in range (2 , n + 1 ):
157
156
if i != k :
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 ))
157
+ model .addCons (quicksum (f [j , i , k ] for j in range (1 , n + 1 ) if (j , i , k ) in f ) == \
158
+ quicksum (f [i , j , k ] for j in range (1 , n + 1 ) if (i , j , k ) in f ),
159
+ "FlowCons(%s,%s)" % (i , k ))
161
160
162
- for (i ,j , k ) in f :
163
- model .addCons (f [i ,j , k ] <= x [i ,j ], "FlowUB(%s,%s,%s)" % (i ,j , k ))
161
+ for (i , j , k ) in f :
162
+ model .addCons (f [i , j , k ] <= x [i , j ], "FlowUB(%s,%s,%s)" % (i , j , k ))
164
163
165
- model .setObjective (quicksum (c [i ,j ] * x [i ,j ] for (i ,j ) in x ), "minimize" )
164
+ model .setObjective (quicksum (c [i , j ] * x [i , j ] for (i , j ) in x ), "minimize" )
166
165
167
- model .data = x ,f
166
+ model .data = x , f
168
167
return model
169
168
170
169
171
-
172
170
def sequence (arcs ):
173
171
"""sequence: make a list of cities to visit, from set of arcs"""
174
172
succ = {}
175
- for (i ,j ) in arcs :
173
+ for (i , j ) in arcs :
176
174
succ [i ] = j
177
- curr = 1 # first node being visited
175
+ curr = 1 # first node being visited
178
176
sol = [curr ]
179
- for i in range (len (arcs )- 2 ):
177
+ for i in range (len (arcs ) - 2 ):
180
178
curr = succ [curr ]
181
179
sol .append (curr )
182
180
return sol
183
181
184
182
185
183
if __name__ == "__main__" :
186
184
n = 5
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 ,
185
+ c = {(1 , 1 ): 0 , (1 , 2 ): 1989 , (1 , 3 ): 102 , (1 , 4 ): 102 , (1 , 5 ): 103 ,
186
+ (2 , 1 ): 104 , (2 , 2 ): 0 , (2 , 3 ): 11 , (2 , 4 ): 104 , (2 , 5 ): 108 ,
187
+ (3 , 1 ): 107 , (3 , 2 ): 108 , (3 , 3 ): 0 , (3 , 4 ): 19 , (3 , 5 ): 102 ,
188
+ (4 , 1 ): 109 , (4 , 2 ): 102 , (4 , 3 ): 107 , (4 , 4 ): 0 , (4 , 5 ): 15 ,
189
+ (5 , 1 ): 13 , (5 , 2 ): 103 , (5 , 3 ): 104 , (5 , 4 ): 101 , (5 , 5 ): 0 ,
192
190
}
193
191
194
- model = mtz (n ,c )
195
- model .hideOutput () # silent mode
192
+ model = mtz (n , c )
193
+ model .hideOutput () # silent mode
196
194
model .optimize ()
197
195
cost = model .getObjVal ()
198
196
print ()
199
197
print ("Miller-Tucker-Zemlin's model:" )
200
198
print ("Optimal value:" , cost )
201
- #model.printAttr("X")
199
+ # model.printAttr("X")
202
200
for v in model .getVars ():
203
201
if model .getVal (v ) > 0.001 :
204
202
print (v .name , "=" , model .getVal (v ))
205
203
206
- x ,u = model .data
207
- sol = [i for (p ,i ) in sorted ([(int (model .getVal (u [i ])+ .5 ),i ) for i in range (1 ,n + 1 )])]
204
+ x , u = model .data
205
+ sol = [i for (p , i ) in sorted ([(int (model .getVal (u [i ]) + .5 ), i ) for i in range (1 , n + 1 )])]
208
206
print (sol )
209
- arcs = [(i ,j ) for (i ,j ) in x if model .getVal (x [i ,j ]) > .5 ]
207
+ arcs = [(i , j ) for (i , j ) in x if model .getVal (x [i , j ]) > .5 ]
210
208
sol = sequence (arcs )
211
209
print (sol )
212
210
# assert cost == 5
213
211
214
- model = mtz_strong (n ,c )
215
- model .hideOutput () # silent mode
212
+ model = mtz_strong (n , c )
213
+ model .hideOutput () # silent mode
216
214
model .optimize ()
217
215
cost = model .getObjVal ()
218
216
print ()
219
217
print ("Miller-Tucker-Zemlin's model with stronger constraints:" )
220
- print ("Optimal value:" ,cost )
221
- #model.printAttr("X")
218
+ print ("Optimal value:" , cost )
219
+ # model.printAttr("X")
222
220
for v in model .getVars ():
223
221
if model .getVal (v ) > 0.001 :
224
222
print (v .name , "=" , model .getVal (v ))
225
223
226
- x ,u = model .data
227
- sol = [i for (p ,i ) in sorted ([(int (model .getVal (u [i ])+ .5 ),i ) for i in range (1 ,n + 1 )])]
224
+ x , u = model .data
225
+ sol = [i for (p , i ) in sorted ([(int (model .getVal (u [i ]) + .5 ), i ) for i in range (1 , n + 1 )])]
228
226
print (sol )
229
- arcs = [(i ,j ) for (i ,j ) in x if model .getVal (x [i ,j ]) > .5 ]
227
+ arcs = [(i , j ) for (i , j ) in x if model .getVal (x [i , j ]) > .5 ]
230
228
sol = sequence (arcs )
231
229
print (sol )
232
230
# assert cost == 5
233
231
234
- model = scf (n ,c )
235
- model .hideOutput () # silent mode
232
+ model = scf (n , c )
233
+ model .hideOutput () # silent mode
236
234
model .optimize ()
237
235
cost = model .getObjVal ()
238
236
print ()
239
237
print ("Single-commodity flow formulation:" )
240
- print ("Optimal value:" ,cost )
241
- #model.printAttr("X")
238
+ print ("Optimal value:" , cost )
239
+ # model.printAttr("X")
242
240
for v in model .getVars ():
243
241
if model .getVal (v ) > 0.001 :
244
242
print (v .name , "=" , model .getVal (v ))
245
243
246
- x ,f = model .data
247
- arcs = [(i ,j ) for (i ,j ) in x if model .getVal (x [i ,j ]) > .5 ]
244
+ x , f = model .data
245
+ arcs = [(i , j ) for (i , j ) in x if model .getVal (x [i , j ]) > .5 ]
248
246
sol = sequence (arcs )
249
247
print (sol )
250
248
# assert cost == 5
251
249
252
- model = mcf (n ,c )
253
- model .hideOutput () # silent mode
250
+ model = mcf (n , c )
251
+ model .hideOutput () # silent mode
254
252
model .optimize ()
255
253
cost = model .getObjVal ()
256
254
print ()
257
255
print ("Multi-commodity flow formulation:" )
258
- print ("Optimal value:" ,cost )
259
- #model.printAttr("X")
256
+ print ("Optimal value:" , cost )
257
+ # model.printAttr("X")
260
258
for v in model .getVars ():
261
- if model .getVal (v )> 0.001 :
259
+ if model .getVal (v ) > 0.001 :
262
260
print (v .name , "=" , model .getVal (v ))
263
261
264
- x ,f = model .data
265
- arcs = [(i ,j ) for (i ,j ) in x if model .getVal (x [i ,j ]) > .5 ]
262
+ x , f = model .data
263
+ arcs = [(i , j ) for (i , j ) in x if model .getVal (x [i , j ]) > .5 ]
266
264
sol = sequence (arcs )
267
265
print (sol )
268
266
# assert cost == 5
0 commit comments