tsp.py
Go to the documentation of this file.
13
14from __future__ import print_function
15from gams import *
16from itertools import product
17import os
18import sys
19
21 return '''
22\$Title Traveling Salesman Problem Instance with Python
23
24\$Ontext
25
26The sub_tour elimination constraints are generated by a Python
27script. The MIP is solved over and over, but GAMS have to
28generate the model only after n cuts have been added.
29
30\$Offtext
31
32\$if not set tspdata \$abort 'tspdata not set'
33
34set ii cities
35 i(ii) subset of cities
36alias (ii,jj),(i,j,k);
37
38parameter c(ii,jj) distance matrix;
39
40\$gdxin %tspdata%
42
43\$if not set nrCities \$set nrCities 20
44i(ii)\$(ord(ii) < %nrCities%) = yes;
45
46variables x(ii,jj) decision variables - leg of trip
47 z objective variable;
48binary variable x; x.fx(ii,ii) = 0;
49
50equations objective total cost
51 rowsum(ii) leave each city only once
52 colsum(jj) arrive at each city only once;
53
54* the assignment problem is a relaxation of the TSP
55objective.. z =e= sum((i,j), c(i,j)*x(i,j));
56rowsum(i).. sum(j, x(i,j)) =e= 1;
57colsum(j).. sum(i, x(i,j)) =e= 1;
58
59\$if not set cmax \$set cmax 2
60set cut /c0*c%cmax%/;
61parameter
62 acut(cut,ii,jj) cut constraint matrix
63 rhscut(cut) cut constraint rhs;
64
65equation sscut(cut) sub_tour elimination cut;
66sscut(cut).. sum((i,j), Acut(cut,i,j)*x(i,j)) =l= RHScut(cut);
67
68set cc(cut) previous cuts; cc(cut) = no;
69\$if set cutdata execute_load '%cutdata%', cc, Acut, RHScut;
70
71Acut(cut,i,j)\$(not cc(cut)) = eps;
72RHScut(cut)\$(not cc(cut)) = card(ii);
73
74model assign /all/;
75
76option optcr=0; '''
77
78
79if __name__ == "__main__":
80 if len(sys.argv) > 1:
81 ws = GamsWorkspace(system_directory = sys.argv[1])
82 else:
83 ws = GamsWorkspace()
84
85 # number of cuts that can be added to a model instance
86 cuts_per_round = 10
87 # current cut
88 curcut = 0
89 # cut limit for current model instance (cmax = curcut + cuts_per_round)
90 cmax = 0
91
92 # database used to collect all generated cuts
93 cut_data = ws.add_database()
94 cc = cut_data.add_set("cc", 1, "")
95 acut = cut_data.add_parameter("acut", 3, "")
96 rhscut = cut_data.add_parameter("rhscut", 1, "")
97 # list of cities (i1, i2, i3, ...)
98 n = []
99
100 # do-while loop
101 while True:
102 # create a new model instance when the cut limit is reached
103 if curcut >= cmax:
104 cmax = curcut + cuts_per_round
105 cut_data.export();
106
107 # create the checkpoint
108 tsp_job = ws.add_job_from_string(get_model_text())
109 cp = ws.add_checkpoint()
110 opt = ws.add_options()
111 opt.defines["nrcities"] = "20"
112 opt.defines["cmax"] = str(cmax - 1)
113 opt.defines["cutdata"] = cut_data.name
114 opt.defines["tspdata"] = '"' + os.path.join(os.path.dirname(os.path.realpath(__file__)), "..", "Data", "tsp.gdx") + '"'
115 tsp_job.run(gams_options=opt, checkpoint=cp)
116
117 # fill the n list only once
118 if not n:
119 for i in tsp_job.out_db["i"]:
120 n += i.keys
121
122 #instantiate the model instance with modifiers mi_acut and mi_rhscut
123 mi = cp.add_modelinstance()
124 mi_acut = mi.sync_db.add_parameter("acut", 3, "")
125 mi_rhscut = mi.sync_db.add_parameter("rhscut", 1, "")
126 mi.instantiate("assign use mip min z", [GamsModifier(mi_acut), GamsModifier(mi_rhscut)])
127
128 #solve model instance using update type accumulate and clear acut and rhscut afterwards
129 mi.solve(SymbolUpdateType.Accumulate)
130 mi.sync_db["acut"].clear()
131 mi.sync_db["rhscut"].clear()
132
133 # collect graph information from the solution
134 graph = {}
135 not_visited = list(n)
136 for i,j in product(n,n):
137 if mi.sync_db["x"].find_record([i, j]).level > 0.5:
138 graph[i] = j
139
140 # find all subtours and add the required cuts by modifying acut and rhscut
141 while not_visited:
142 i = not_visited[0]
143 sub_tour = [i]
144 while graph[i] != not_visited[0]:
145 i = graph[i]
146 sub_tour.append(i)
147 not_visited = list(filter(lambda x: x not in sub_tour, not_visited))
148
149 # add the cuts to both databases (cut_data, mi.sync_db)
150 for i,j in product(sub_tour,sub_tour):
151 acut.add_record(["c"+str(curcut), i, j]).value = 1
152 mi_acut.add_record(["c"+str(curcut), i, j]).value = 1
153 rhscut.add_record("c"+str(curcut)).value = len(sub_tour)-0.5
154 mi_rhscut.add_record(["c"+str(curcut)]).value = len(sub_tour)-0.5