Loading...
Searching...
No Matches
tsp.py
Go to the documentation of this file.
12
13from itertools import product
14import os
15import sys
16from gams import GamsWorkspace, GamsModifier, SymbolUpdateType
17
18GAMS_MODEL = """
19$if not set tspdata $abort 'tspdata not set'
20
21Set
22 ii 'cities'
23 i(ii) 'subset of cities';
24
25Alias (ii,jj), (i,j,k);
26
27Parameter c(ii,jj) 'distance matrix';
28
29$gdxIn %tspdata%
30$load ii c
31$gdxIn
32
33$if not set nrCities $set nrCities 20
34i(ii)$(ord(ii) < %nrCities%) = yes;
35
36Variable
37 x(ii,jj) 'decision variables - leg of trip'
38 z 'objective variable';
39
40Binary Variable x;
41x.fx(ii,ii) = 0;
42
43Equation
44 objective 'total cost'
45 rowsum(ii) 'leave each city only once'
46 colsum(jj) 'arrive at each city only once';
47
48* the assignment problem is a relaxation of the TSP
49objective.. z =e= sum((i,j), c(i,j)*x(i,j));
50
51rowsum(i).. sum(j, x(i,j)) =e= 1;
52
53colsum(j).. sum(i, x(i,j)) =e= 1;
54
55$if not set cmax $set cmax 2
56set cut /c0*c%cmax%/;
57
58Parameter
59 acut(cut,ii,jj) 'cut constraint matrix'
60 rhscut(cut) 'cut constraint rhs';
61
62Equation sscut(cut) 'sub_tour elimination cut';
63sscut(cut).. sum((i,j), acut(cut,i,j)*x(i,j)) =l= rhscut(cut);
64
65set cc(cut) 'previous cuts';
66cc(cut) = no;
67
68$if set cutdata execute_load '%cutdata%', cc, acut, rhscut;
69
70acut(cut,i,j)$(not cc(cut)) = eps;
71
72rhscut(cut)$(not cc(cut)) = card(ii);
73
74Model assign /all/;
75
76option optcr=0;
77"""
78
79if __name__ == "__main__":
80 sys_dir = sys.argv[1] if len(sys.argv) > 1 else None
81 ws = GamsWorkspace(system_directory=sys_dir)
82
83 # number of cuts that can be added to a model instance
84 cuts_per_round = 10
85
86 # current cut
87 cur_cut = 0
88
89 # cut limit for current model instance (cmax = cur_cut + 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
98 n = [] # list of cities (i1, i2, i3, ...)
99 while True:
100 # create a new model instance if the cut limit is reached
101 if cur_cut >= cmax:
102 cmax = cur_cut + cuts_per_round
103 cut_data.export()
104
105 # create the checkpoint
106 job = ws.add_job_from_string(GAMS_MODEL)
107 cp = ws.add_checkpoint()
108 opt = ws.add_options()
109 opt.defines["nrcities"] = "20"
110 opt.defines["cmax"] = str(cmax - 1)
111 opt.defines["cutdata"] = cut_data.name
112 opt.defines["tspdata"] = (
113 '"'
114 + os.path.join(ws.system_directory, "apifiles", "Data", "tsp.gdx")
115 + '"'
116 )
117 job.run(gams_options=opt, checkpoint=cp)
118
119 # fill the cities list only once
120 if not n:
121 for i in job.out_db["i"]:
122 n += i.keys
123
124 # instantiate the model instance with modifiers mi_acut and mi_rhscut
125 mi = cp.add_modelinstance()
126 mi_acut = mi.sync_db.add_parameter("acut", 3)
127 mi_rhscut = mi.sync_db.add_parameter("rhscut", 1)
128 mi.instantiate(
129 "assign use mip min z", [GamsModifier(mi_acut), GamsModifier(mi_rhscut)]
130 )
131
132 # solve model instance using update type accumulate and clear acut and rhscut afterwards
133 mi.solve(SymbolUpdateType.Accumulate)
134 mi.sync_db["acut"].clear()
135 mi.sync_db["rhscut"].clear()
136
137 # collect graph information from the solution
138 graph = {}
139 for i, j in product(n, n):
140 if mi.sync_db["x"].find_record([i, j]).level > 0.5:
141 graph[i] = j
142
143 # find all sub tours and add the required cuts by modifying acut and rhscut
144 not_visited = list(n) # make copy
145 while not_visited:
146 i = not_visited[0]
147 sub_tour = [i]
148 while graph[i] != not_visited[0]:
149 i = graph[i]
150 sub_tour.append(i)
151 not_visited = list(filter(lambda x: x not in sub_tour, not_visited))
152
153 # add the cuts to both databases (cut_data, mi.sync_db)
154 for i, j in product(sub_tour, sub_tour):
155 acut.add_record([f"c{cur_cut}", i, j]).value = 1
156 mi_acut.add_record([f"c{cur_cut}", i, j]).value = 1
157 rhscut.add_record(f"c{cur_cut}").value = len(sub_tour) - 0.5
158 mi_rhscut.add_record([f"c{cur_cut}"]).value = len(sub_tour) - 0.5
159 cc.add_record(f"c{cur_cut}")
160 cur_cut += 1
161
162 # quit the loop if a sub tour contains all cities and print the solution
163 if len(sub_tour) == len(n):
164 print(f"z={mi.sync_db['z'].first_record().level}")
165 print("sub_tour: ")
166 for i in sub_tour:
167 print(f"{i} -> ", end="")
168 print(sub_tour[0])
169 break
GamsWorkspace rhscut
Definition: tsp.py:96
GamsWorkspace cc
Definition: tsp.py:94
list i
Definition: tsp.py:146
GamsWorkspace acut
Definition: tsp.py:95