11##@file atsp.py
2- #@brief solve the asymmetric traveling salesman problem
2+ # @brief solve the asymmetric traveling salesman problem
33
44"""
55
1111
1212Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
1313"""
14- from pyscipopt import Model , quicksum , multidict
14+ from pyscipopt import Model , quicksum
1515
16- def mtz (n ,c ):
16+
17+ def mtz (n , c ):
1718 """mtz: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
1819 (potential formulation)
1920 Parameters:
@@ -24,30 +25,29 @@ def mtz(n,c):
2425
2526 model = Model ("atsp - mtz" )
2627
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 ):
3132 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 ))
3334
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 )
3738
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 ):
4041 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 ))
4243
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" )
4745
46+ model .data = x , u
47+ return model
4848
4949
50- def mtz_strong (n ,c ):
50+ def mtz_strong (n , c ):
5151 """mtz_strong: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
5252 (potential formulation, adding stronger constraints)
5353 Parameters:
@@ -57,34 +57,34 @@ def mtz_strong(n,c):
5757 """
5858
5959 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 ):
6565 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 ))
6767
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 )
7171
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 ):
7474 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 )
7680
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" )
8082
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
8484 return model
8585
8686
87- def scf (n ,c ):
87+ def scf (n , c ):
8888 """scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem
8989 Parameters:
9090 - n: number of nodes
@@ -93,40 +93,39 @@ def scf(n,c):
9393 """
9494 model = Model ("atsp - scf" )
9595
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 ):
9999 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 ))
103103 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 ))
105105
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 )
109109
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" )
111111
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 )
115115
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 ):
119119 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 ))
121121
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" )
123123
124- model .data = x ,f
124+ model .data = x , f
125125 return model
126126
127127
128-
129- def mcf (n ,c ):
128+ def mcf (n , c ):
130129 """mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem
131130 Parameters:
132131 - n: number of nodes
@@ -135,134 +134,133 @@ def mcf(n,c):
135134 """
136135 model = Model ("mcf" )
137136
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 ):
141140 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 ))
143142 if i != j and j != 1 :
144- for k in range (2 ,n + 1 ):
143+ for k in range (2 , n + 1 ):
145144 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 ))
147146
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 )
151150
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 )
155154
156- for i in range (2 ,n + 1 ):
155+ for i in range (2 , n + 1 ):
157156 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 ))
161160
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 ))
164163
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" )
166165
167- model .data = x ,f
166+ model .data = x , f
168167 return model
169168
170169
171-
172170def sequence (arcs ):
173171 """sequence: make a list of cities to visit, from set of arcs"""
174172 succ = {}
175- for (i ,j ) in arcs :
173+ for (i , j ) in arcs :
176174 succ [i ] = j
177- curr = 1 # first node being visited
175+ curr = 1 # first node being visited
178176 sol = [curr ]
179- for i in range (len (arcs )- 2 ):
177+ for i in range (len (arcs ) - 2 ):
180178 curr = succ [curr ]
181179 sol .append (curr )
182180 return sol
183181
184182
185183if __name__ == "__main__" :
186184 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 ,
192190 }
193191
194- model = mtz (n ,c )
195- model .hideOutput () # silent mode
192+ model = mtz (n , c )
193+ model .hideOutput () # silent mode
196194 model .optimize ()
197195 cost = model .getObjVal ()
198196 print ()
199197 print ("Miller-Tucker-Zemlin's model:" )
200198 print ("Optimal value:" , cost )
201- #model.printAttr("X")
199+ # model.printAttr("X")
202200 for v in model .getVars ():
203201 if model .getVal (v ) > 0.001 :
204202 print (v .name , "=" , model .getVal (v ))
205203
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 )])]
208206 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 ]
210208 sol = sequence (arcs )
211209 print (sol )
212210 # assert cost == 5
213211
214- model = mtz_strong (n ,c )
215- model .hideOutput () # silent mode
212+ model = mtz_strong (n , c )
213+ model .hideOutput () # silent mode
216214 model .optimize ()
217215 cost = model .getObjVal ()
218216 print ()
219217 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")
222220 for v in model .getVars ():
223221 if model .getVal (v ) > 0.001 :
224222 print (v .name , "=" , model .getVal (v ))
225223
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 )])]
228226 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 ]
230228 sol = sequence (arcs )
231229 print (sol )
232230 # assert cost == 5
233231
234- model = scf (n ,c )
235- model .hideOutput () # silent mode
232+ model = scf (n , c )
233+ model .hideOutput () # silent mode
236234 model .optimize ()
237235 cost = model .getObjVal ()
238236 print ()
239237 print ("Single-commodity flow formulation:" )
240- print ("Optimal value:" ,cost )
241- #model.printAttr("X")
238+ print ("Optimal value:" , cost )
239+ # model.printAttr("X")
242240 for v in model .getVars ():
243241 if model .getVal (v ) > 0.001 :
244242 print (v .name , "=" , model .getVal (v ))
245243
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 ]
248246 sol = sequence (arcs )
249247 print (sol )
250248 # assert cost == 5
251249
252- model = mcf (n ,c )
253- model .hideOutput () # silent mode
250+ model = mcf (n , c )
251+ model .hideOutput () # silent mode
254252 model .optimize ()
255253 cost = model .getObjVal ()
256254 print ()
257255 print ("Multi-commodity flow formulation:" )
258- print ("Optimal value:" ,cost )
259- #model.printAttr("X")
256+ print ("Optimal value:" , cost )
257+ # model.printAttr("X")
260258 for v in model .getVars ():
261- if model .getVal (v )> 0.001 :
259+ if model .getVal (v ) > 0.001 :
262260 print (v .name , "=" , model .getVal (v ))
263261
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 ]
266264 sol = sequence (arcs )
267265 print (sol )
268266 # assert cost == 5
0 commit comments