PySCIPOpt  4.3.0
Python Interface for the SCIP Optimization Suite
atsp.py
Go to the documentation of this file.
1 
3 
4 """
5 
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
11 
12 Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
13 """
14 from pyscipopt import Model, quicksum, multidict
15 
16 def mtz(n,c):
17  """mtz: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
18  (potential formulation)
19  Parameters:
20  - n: number of nodes
21  - c[i,j]: cost for traversing arc (i,j)
22  Returns a model, ready to be solved.
23  """
24 
25  model = Model("atsp - mtz")
26 
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):
31  if i != j:
32  x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
33 
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)
37 
38  for i in range(1,n+1):
39  for j in range(2,n+1):
40  if i != j:
41  model.addCons(u[i] - u[j] + (n-1)*x[i,j] <= n-2, "MTZ(%s,%s)"%(i,j))
42 
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
47 
48 
49 
50 def mtz_strong(n,c):
51  """mtz_strong: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
52  (potential formulation, adding stronger constraints)
53  Parameters:
54  n - number of nodes
55  c[i,j] - cost for traversing arc (i,j)
56  Returns a model, ready to be solved.
57  """
58 
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):
65  if i != j:
66  x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
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)
71 
72  for i in range(1,n+1):
73  for j in range(2,n+1):
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))
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)
80 
81  model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")
82 
83  model.data = x,u
84  return model
85 
86 
87 def scf(n,c):
88  """scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem
89  Parameters:
90  - n: number of nodes
91  - c[i,j]: cost for traversing arc (i,j)
92  Returns a model, ready to be solved.
93  """
94  model = Model("atsp - scf")
95 
96  x,f = {},{}
97  for i in range(1,n+1):
98  for j in range(1,n+1):
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))
103  else:
104  f[i,j] = model.addVar(lb=0, ub=n-2, vtype="C", name="f(%s,%s)"%(i,j))
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)
109 
110  model.addCons(quicksum(f[1,j] for j in range(2,n+1)) == n-1, "FlowOut")
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)
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):
119  if i != j:
120  model.addCons(f[i,j] <= (n-2)*x[i,j], "FlowUB(%s,%s)"%(i,j))
121 
122  model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")
123 
124  model.data = x,f
125  return model
126 
127 
128 
129 def mcf(n,c):
130  """mcf: multi-commodity flow formulation for the (asymmetric) traveling salesman problem
131  Parameters:
132  - n: number of nodes
133  - c[i,j]: cost for traversing arc (i,j)
134  Returns a model, ready to be solved.
135  """
136  model = Model("mcf")
137 
138  x,f = {},{}
139  for i in range(1,n+1):
140  for j in range(1,n+1):
141  if i != j:
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):
145  if i != k:
146  f[i,j,k] = model.addVar(ub=1, vtype="C", name="f(%s,%s,%s)"%(i,j,k))
147 
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)
151 
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)
155 
156  for i in range(2,n+1):
157  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))
161 
162  for (i,j,k) in f:
163  model.addCons(f[i,j,k] <= x[i,j], "FlowUB(%s,%s,%s)"%(i,j,k))
164 
165  model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")
166 
167  model.data = x,f
168  return model
169 
170 
171 
172 def sequence(arcs):
173  """sequence: make a list of cities to visit, from set of arcs"""
174  succ = {}
175  for (i,j) in arcs:
176  succ[i] = j
177  curr = 1 # first node being visited
178  sol = [curr]
179  for i in range(len(arcs)-2):
180  curr = succ[curr]
181  sol.append(curr)
182  return sol
183 
184 
185 if __name__ == "__main__":
186  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,
192  }
193 
194  model = mtz(n,c)
195  model.hideOutput() # silent mode
196  model.optimize()
197  cost = model.getObjVal()
198  print()
199  print("Miller-Tucker-Zemlin's model:")
200  print("Optimal value:", cost)
201  #model.printAttr("X")
202  for v in model.getVars():
203  if model.getVal(v) > 0.001:
204  print(v.name, "=", model.getVal(v))
205 
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)])]
208  print(sol)
209  arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5]
210  sol = sequence(arcs)
211  print(sol)
212  # assert cost == 5
213 
214  model = mtz_strong(n,c)
215  model.hideOutput() # silent mode
216  model.optimize()
217  cost = model.getObjVal()
218  print()
219  print("Miller-Tucker-Zemlin's model with stronger constraints:")
220  print("Optimal value:",cost)
221  #model.printAttr("X")
222  for v in model.getVars():
223  if model.getVal(v) > 0.001:
224  print(v.name, "=", model.getVal(v))
225 
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)])]
228  print(sol)
229  arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5]
230  sol = sequence(arcs)
231  print(sol)
232  # assert cost == 5
233 
234  model = scf(n,c)
235  model.hideOutput() # silent mode
236  model.optimize()
237  cost = model.getObjVal()
238  print()
239  print("Single-commodity flow formulation:")
240  print("Optimal value:",cost)
241  #model.printAttr("X")
242  for v in model.getVars():
243  if model.getVal(v) > 0.001:
244  print(v.name, "=", model.getVal(v))
245 
246  x,f = model.data
247  arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5]
248  sol = sequence(arcs)
249  print(sol)
250  # assert cost == 5
251 
252  model = mcf(n,c)
253  model.hideOutput() # silent mode
254  model.optimize()
255  cost = model.getObjVal()
256  print()
257  print("Multi-commodity flow formulation:")
258  print("Optimal value:",cost)
259  #model.printAttr("X")
260  for v in model.getVars():
261  if model.getVal(v)>0.001:
262  print(v.name, "=", model.getVal(v))
263 
264  x,f = model.data
265  arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5]
266  sol = sequence(arcs)
267  print(sol)
268  # assert cost == 5
atsp.mtz_strong
def mtz_strong(n, c)
Definition: atsp.py:50
pyscipopt.expr.quicksum
def quicksum(termlist)
Definition: expr.pxi:357
atsp.mtz
def mtz(n, c)
Definition: atsp.py:16
atsp.sequence
def sequence(arcs)
Definition: atsp.py:172
atsp.scf
def scf(n, c)
Definition: atsp.py:87
atsp.mcf
def mcf(n, c)
Definition: atsp.py:129