Optimization: principles and algorithms, by Michel Bierlaire
tspExact.m
Go to the documentation of this file.
1 %> \file
2 %> Write the Traveling Salesman Problem (TSP) as an integer optimization problem in standard form. See Section 25.2.3 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref run2512.m
5 %> @note Tested with \ref run2703.m
6 %>
7 %> @author Michel Bierlaire
8 %> @date Tue Apr 14 10:42:27 2015
9 %> @ingroup Algorithms
10 %> @ingroup chap25
11 
12 %> Write the TSP as an integer optimization problem in standard form
13 %> @param dist distance matrix
14 %> @return [A,b,c]
15 
16 function tour = tspExact(dist)
17  n = rows(dist) ;
18  if (columns(dist) != n)
19  error("Distance matrix must be square") ;
20  endif
21 
22 % Numbering of variables
23 % 1--> n(n-1) x_ij
24  vtype1 = repmat('I',1,n*(n-1)) ;
25 % n(n-1) + 1 --> n(n-1) + n-1 u_i (for the subtour elimination constraint)
26  vtype2 = repmat('I',1,n-1) ;
27 % n(n-1) + n-1 + 1 --> n(n-1) + n-1 + n(n-1) slack variables for x_ij <= 1
28  vtype3 = repmat('C',1,n*(n-1)) ;
29 % n(n-1) + n-1 + n(n-1) + 1 --> n(n-1) + n-1 + n(n-1) + (n-1)*(n-2) slack variables for subrout constraints
30  vtype4 = repmat('C',1,(n-1)*(n-2)) ;
31  vtype = [vtype1 vtype2 vtype3 vtype4] ;
32 
33  indices = zeros(size(dist)) ;
34  index = 1 ;
35  for i = 1:n
36  for j = 1:n
37  if i != j
38  indices(i,j) = index ;
39  index += 1 ;
40  endif
41  endfor
42  endfor
43  numberOfVariables = (n-1) * (3*n-1) ;
44 
45 % Objective function
46 
47  c = zeros(numberOfVariables,1) ;
48  for i = 1:n
49  for j = 1:n
50  if i != j
51  c(indices(i,j))= dist(i,j) ;
52  endif
53  endfor
54  endfor
55 
56 
57 % Each city has exactly one predecessor
58  A1 = zeros(n,numberOfVariables) ;
59  b1 = ones(n,1) ;
60  ctype1 = repmat('S',1,n) ;
61  for j = 1:n
62  for i = 1:n
63  if i != j
64  A1(j,indices(i,j)) = 1 ;
65  endif
66  endfor
67  endfor
68 
69 % Each city has exactly one successor
70  A2 = zeros(n,numberOfVariables) ;
71  b2 = ones(n,1) ;
72  ctype2 = repmat('S',1,n) ;
73  for i = 1:n
74  for j = 1:n
75  if i != j
76  A2(i,indices(i,j)) = 1 ;
77  endif
78  endfor
79  endfor
80 
81 % Subtour elimination
82 
83  A3 = zeros((n-1)*(n-2),numberOfVariables) ;
84  b3 = zeros((n-1)*(n-2),1) ;
85  ctype3 = repmat('S',1,(n-1)*(n-2)) ;
86  index = 1 ;
87  for i = 2:n
88  for j = 2:n
89  if i != j
90  A3(index,n*(n-1) + i-1) = 1 ;
91  A3(index,n*(n-1) + j-1) = -1 ;
92  A3(index,indices(i,j)) = n-1 ;
93  A3(index,n*(n-1) + n-1 + n*(n-1) + index) = 1 ;
94  b3(index) = n-2 ;
95  index += 1 ;
96  endif
97  endfor
98  endfor
99 
100 % x_ij <= 1
101 
102  A4 = zeros(n*(n-1),numberOfVariables) ;
103  b4 = ones(n*(n-1),1) ;
104  ctype4 = repmat('S',1,n*(n-1)) ;
105  for i = 1:n
106  for j = 1:n
107  if i != j
108  A4(indices(i,j),indices(i,j)) = 1 ;
109  A4(indices(i,j),n * (n-1) + n + indices(i,j)-1) = 1 ;
110  endif
111  endfor
112  endfor
113 
114 
115  A = [A1 ; A2 ; A3 ; A4] ;
116  b = [b1 ; b2 ; b3 ; b4] ;
117 
118  [mm,nn] = size(A) ;
119  lb = zeros(numberOfVariables,1) ;
120  ub = 1000*ones(numberOfVariables,1) ;
121  ctype = [ctype1 ctype2 ctype3 ctype4] ;
122  [x_min,tot_cost,errnum,extra] = glpk(c,A,b,lb,ub,ctype,vtype,1) ;
123  for i=1:n
124  for j=1:n
125  if (i != j)
126  if (x_min(indices(i,j)) == 1)
127  printf("%d --> %d\n",i,j)
128  endif
129  endif
130  endfor
131  endfor
132 
133  printf("%d cities\n",n) ;
134  printf("Total cost: %f\n",tot_cost) ;
135  tour = [1] ;
136  currentcity = 1 ;
137  do
138  for j = 1:n
139  if currentcity != j
140  if (x_min(indices(currentcity,j)) == 1)
141  currentcity = j ;
142  tour = [tour ; currentcity ] ;
143  break ;
144  endif
145  endif
146  endfor
147  until(currentcity == 1)
148 
149 
150 endfunction
151 
tspExact
function tspExact(in dist)
Write the TSP as an integer optimization problem in standard form.
Copyright 2015-2018 Michel Bierlaire