611 巡回セールスマン問題¶
Inspired by http://www.gurobi.com/documentation/8.0/examples/tsp_py.html
読み込み¶
In [1]:
import pandas as pd
from geopy.distance import great_circle
import gmaps
google_map_api_key = 'hogehoge'
gmaps.configure(api_key=google_map_api_key)
from ipywidgets.embed import embed_minimal_html
import gurobipy as gp
関数定義¶
In [2]:
def calc_dist_from_latlng(lat1, lng1, lat2, lng2, unit='m'):
if unit == 'm':
return great_circle((lat1, lng1), (lat2, lng2)).meters
elif unit == 'km':
return great_circle((lat1, lng1), (lat2, lng2)).kilometers
else:
print('ERROR: unit')
return False
パラメータ定義¶
Latlag定義¶
In [3]:
IND_BUILDINGS = ['Office', 'Library', 'Class', 'Gym', 'Pool', 'Ramen']
latlng_df = pd.DataFrame([[35.706048, 139.706735],
[35.705613, 139.707174],
[35.705713, 139.705254],
[35.706034, 139.716454],
[35.704950, 139.708149],
[35.708074, 139.708828]],
index=IND_BUILDINGS, columns=['lat', 'lng'])
dist_df = pd.DataFrame(index=IND_BUILDINGS, columns=IND_BUILDINGS)
for key, values in dist_df.iterrows():
for key2, value in values.items():
if key == key2:
#break
pass
dist_df.loc[key, key2] = calc_dist_from_latlng(latlng_df.loc[key].lat, latlng_df.loc[key].lng,
latlng_df.loc[key2].lat, latlng_df.loc[key2].lng)
display(latlng_df)
print('Distance matrix')
display(dist_df)
| lat | lng | |
|---|---|---|
| Office | 35.706048 | 139.706735 |
| Library | 35.705613 | 139.707174 |
| Class | 35.705713 | 139.705254 |
| Gym | 35.706034 | 139.716454 |
| Pool | 35.704950 | 139.708149 |
| Ramen | 35.708074 | 139.708828 |
Distance matrix
| Office | Library | Class | Gym | Pool | Ramen | |
|---|---|---|---|---|---|---|
| Office | 0 | 62.5369 | 138.815 | 877.558 | 176.656 | 294.05 |
| Library | 62.5369 | 0 | 173.719 | 839.227 | 114.828 | 311.75 |
| Class | 138.815 | 173.719 | 0 | 1011.91 | 274.824 | 416.005 |
| Gym | 877.558 | 839.227 | 1011.91 | 0 | 759.513 | 724.967 |
| Pool | 176.656 | 114.828 | 274.824 | 759.513 | 0 | 352.742 |
| Ramen | 294.05 | 311.75 | 416.005 | 724.967 | 352.742 | 0 |
Mixed Integer Linear programming of Travelling salesman problem¶
TSP with MTZ[1]
\begin{align}
\text{Minimize } \Sigma_{i=1}^n\Sigma_{j=1}^n c_{ij}\delta_{ij}, \\
\text{Subject to }\Sigma_{i=1,i\neq j}^n\delta_{ij}&=1 &\forall j\in \{1,...,n\}, \\
\Sigma_{j=1,j\neq i}^n\delta_{ji}&=1 &\forall i\in \{1,...,n\}, \\
x_{i} - x_{j} + (n-1)\delta_{ij} &\le n-2 &\forall i\in \{1,...,n\}, \forall j \in \{2,...,n\} \\
\text{where } c \in R^{n\times n}, \delta\in \{0,1\}^{n\times n}, 1\le x_{i}\le n-1 \nonumber.
\end{align}
See reference for more detailed formulation.
[1] C. E. Miller, A. W. Tucker, and R. A. Zemlin, “Integer Programming Formulation of Traveling Salesman Problems,” J. ACM, vol. 7, no. 4, pp. 326–329, 1960.
[2] A. Langevin, F. Soumis, and J. Desrosiers, “Classification of travelling salesman problem formulations,” Oper. Res. Lett., vol. 9, no. 2, pp. 127–132, 1990.
[3] A. Subramanyam and C. E. Gounaris, “A branch-and-cut framework for the consistent traveling salesman problem,” Eur. J. Oper. Res., vol. 248, no. 2, pp. 384–395, 2016.
In [4]:
model = gp.Model('tsp')
delta_var_dict = {}
cont_var_dict = {}
for src in dist_df.keys():
if src != dist_df.keys()[0]:
cont_var_dict[src] = model.addVar(lb=1, ub=len(dist_df.keys())-1, name='u({})'.format(src))
for dest in dist_df.keys():
if src == dest: continue
delta_var_dict[(src, dest)] = model.addVar(obj=round(dist_df.loc[src, dest]), vtype=gp.GRB.BINARY,
name='b_edge({},{})'.format(src, dest))
model.update()
for src in dist_df.keys():
lhs_sd = gp.LinExpr(0)
lhs_ds = gp.LinExpr(0)
for dest in dist_df.keys():
if src == dest: continue
lhs_sd += delta_var_dict[(src, dest)]
lhs_ds += delta_var_dict[(dest, src)]
model.addConstr(lhs_sd == 1, name='eq_src_dest_({}{})'.format(src, dest))
model.addConstr(lhs_ds == 1, name='eq_dest_src_({}{})'.format(dest, src))
for src in dist_df.keys()[1:]:
for dest in dist_df.keys()[1:]:
if src == dest: continue
model.addConstr(cont_var_dict[src] - cont_var_dict[dest]
+ (len(dist_df.keys())-1)*delta_var_dict[(src, dest)]
<= len(dist_df.keys())-2, name='eq_u({},{})'.format(src, dest))
output_file = 'tsp.lp'
model.write(output_file)
print(open(output_file).read())
\ Model tsp
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
63 b_edge(Office,Library) + 139 b_edge(Office,Class)
+ 878 b_edge(Office,Gym) + 177 b_edge(Office,Pool)
+ 294 b_edge(Office,Ramen) + 63 b_edge(Library,Office)
+ 174 b_edge(Library,Class) + 839 b_edge(Library,Gym)
+ 115 b_edge(Library,Pool) + 312 b_edge(Library,Ramen)
+ 139 b_edge(Class,Office) + 174 b_edge(Class,Library)
+ 1012 b_edge(Class,Gym) + 275 b_edge(Class,Pool)
+ 416 b_edge(Class,Ramen) + 878 b_edge(Gym,Office)
+ 839 b_edge(Gym,Library) + 1012 b_edge(Gym,Class)
+ 760 b_edge(Gym,Pool) + 725 b_edge(Gym,Ramen) + 177 b_edge(Pool,Office)
+ 115 b_edge(Pool,Library) + 275 b_edge(Pool,Class)
+ 760 b_edge(Pool,Gym) + 353 b_edge(Pool,Ramen)
+ 294 b_edge(Ramen,Office) + 312 b_edge(Ramen,Library)
+ 416 b_edge(Ramen,Class) + 725 b_edge(Ramen,Gym)
+ 353 b_edge(Ramen,Pool)
Subject To
eq_src_dest_(OfficeRamen): b_edge(Office,Library) + b_edge(Office,Class)
+ b_edge(Office,Gym) + b_edge(Office,Pool) + b_edge(Office,Ramen) = 1
eq_dest_src_(RamenOffice): b_edge(Library,Office) + b_edge(Class,Office)
+ b_edge(Gym,Office) + b_edge(Pool,Office) + b_edge(Ramen,Office) = 1
eq_src_dest_(LibraryRamen): b_edge(Library,Office) + b_edge(Library,Class)
+ b_edge(Library,Gym) + b_edge(Library,Pool) + b_edge(Library,Ramen)
= 1
eq_dest_src_(RamenLibrary): b_edge(Office,Library) + b_edge(Class,Library)
+ b_edge(Gym,Library) + b_edge(Pool,Library) + b_edge(Ramen,Library)
= 1
eq_src_dest_(ClassRamen): b_edge(Class,Office) + b_edge(Class,Library)
+ b_edge(Class,Gym) + b_edge(Class,Pool) + b_edge(Class,Ramen) = 1
eq_dest_src_(RamenClass): b_edge(Office,Class) + b_edge(Library,Class)
+ b_edge(Gym,Class) + b_edge(Pool,Class) + b_edge(Ramen,Class) = 1
eq_src_dest_(GymRamen): b_edge(Gym,Office) + b_edge(Gym,Library)
+ b_edge(Gym,Class) + b_edge(Gym,Pool) + b_edge(Gym,Ramen) = 1
eq_dest_src_(RamenGym): b_edge(Office,Gym) + b_edge(Library,Gym)
+ b_edge(Class,Gym) + b_edge(Pool,Gym) + b_edge(Ramen,Gym) = 1
eq_src_dest_(PoolRamen): b_edge(Pool,Office) + b_edge(Pool,Library)
+ b_edge(Pool,Class) + b_edge(Pool,Gym) + b_edge(Pool,Ramen) = 1
eq_dest_src_(RamenPool): b_edge(Office,Pool) + b_edge(Library,Pool)
+ b_edge(Class,Pool) + b_edge(Gym,Pool) + b_edge(Ramen,Pool) = 1
eq_src_dest_(RamenRamen): b_edge(Ramen,Office) + b_edge(Ramen,Library)
+ b_edge(Ramen,Class) + b_edge(Ramen,Gym) + b_edge(Ramen,Pool) = 1
eq_dest_src_(RamenRamen): b_edge(Office,Ramen) + b_edge(Library,Ramen)
+ b_edge(Class,Ramen) + b_edge(Gym,Ramen) + b_edge(Pool,Ramen) = 1
eq_u(Library,Class): u(Library) + 5 b_edge(Library,Class) - u(Class) <= 4
eq_u(Library,Gym): u(Library) + 5 b_edge(Library,Gym) - u(Gym) <= 4
eq_u(Library,Pool): u(Library) + 5 b_edge(Library,Pool) - u(Pool) <= 4
eq_u(Library,Ramen): u(Library) + 5 b_edge(Library,Ramen) - u(Ramen) <= 4
eq_u(Class,Library): - u(Library) + u(Class) + 5 b_edge(Class,Library)
<= 4
eq_u(Class,Gym): u(Class) + 5 b_edge(Class,Gym) - u(Gym) <= 4
eq_u(Class,Pool): u(Class) + 5 b_edge(Class,Pool) - u(Pool) <= 4
eq_u(Class,Ramen): u(Class) + 5 b_edge(Class,Ramen) - u(Ramen) <= 4
eq_u(Gym,Library): - u(Library) + u(Gym) + 5 b_edge(Gym,Library) <= 4
eq_u(Gym,Class): - u(Class) + u(Gym) + 5 b_edge(Gym,Class) <= 4
eq_u(Gym,Pool): u(Gym) + 5 b_edge(Gym,Pool) - u(Pool) <= 4
eq_u(Gym,Ramen): u(Gym) + 5 b_edge(Gym,Ramen) - u(Ramen) <= 4
eq_u(Pool,Library): - u(Library) + u(Pool) + 5 b_edge(Pool,Library) <= 4
eq_u(Pool,Class): - u(Class) + u(Pool) + 5 b_edge(Pool,Class) <= 4
eq_u(Pool,Gym): - u(Gym) + u(Pool) + 5 b_edge(Pool,Gym) <= 4
eq_u(Pool,Ramen): u(Pool) + 5 b_edge(Pool,Ramen) - u(Ramen) <= 4
eq_u(Ramen,Library): - u(Library) + u(Ramen) + 5 b_edge(Ramen,Library)
<= 4
eq_u(Ramen,Class): - u(Class) + u(Ramen) + 5 b_edge(Ramen,Class) <= 4
eq_u(Ramen,Gym): - u(Gym) + u(Ramen) + 5 b_edge(Ramen,Gym) <= 4
eq_u(Ramen,Pool): - u(Pool) + u(Ramen) + 5 b_edge(Ramen,Pool) <= 4
Bounds
1 <= u(Library) <= 5
1 <= u(Class) <= 5
1 <= u(Gym) <= 5
1 <= u(Pool) <= 5
1 <= u(Ramen) <= 5
Binaries
b_edge(Office,Library) b_edge(Office,Class) b_edge(Office,Gym)
b_edge(Office,Pool) b_edge(Office,Ramen) b_edge(Library,Office)
b_edge(Library,Class) b_edge(Library,Gym) b_edge(Library,Pool)
b_edge(Library,Ramen) b_edge(Class,Office) b_edge(Class,Library)
b_edge(Class,Gym) b_edge(Class,Pool) b_edge(Class,Ramen)
b_edge(Gym,Office) b_edge(Gym,Library) b_edge(Gym,Class) b_edge(Gym,Pool)
b_edge(Gym,Ramen) b_edge(Pool,Office) b_edge(Pool,Library)
b_edge(Pool,Class) b_edge(Pool,Gym) b_edge(Pool,Ramen)
b_edge(Ramen,Office) b_edge(Ramen,Library) b_edge(Ramen,Class)
b_edge(Ramen,Gym) b_edge(Ramen,Pool)
End
最適化実行¶
In [5]:
model.optimize()
Optimize a model with 32 rows, 35 columns and 120 nonzeros
Variable types: 5 continuous, 30 integer (30 binary)
Coefficient statistics:
Matrix range [1e+00, 5e+00]
Objective range [6e+01, 1e+03]
Bounds range [1e+00, 5e+00]
RHS range [1e+00, 4e+00]
Found heuristic solution: objective 2837.0000000
Presolve time: 0.00s
Presolved: 32 rows, 35 columns, 180 nonzeros
Variable types: 5 continuous, 30 integer (30 binary)
Root relaxation: objective 2.061517e+03, 21 iterations, 0.00 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 2061.51724 0 9 2837.00000 2061.51724 27.3% - 0s
H 0 0 2495.0000000 2061.51724 17.4% - 0s
H 0 0 2387.0000000 2061.51724 13.6% - 0s
H 0 0 2207.0000000 2061.51724 6.59% - 0s
0 0 infeasible 0 2207.00000 2207.00000 0.00% - 0s
Cutting planes:
Gomory: 1
Implied bound: 2
Clique: 2
MIR: 6
Explored 1 nodes (31 simplex iterations) in 0.06 seconds
Thread count was 8 (of 8 available processors)
Solution count 4: 2207 2387 2495 2837
Optimal solution found (tolerance 1.00e-04)
Best objective 2.207000000000e+03, best bound 2.207000000000e+03, gap 0.0000%
最適化結果¶
In [6]:
optimum_var_dict = model.getAttr('x', delta_var_dict)
optimal_dict = {}
for src in dist_df.keys():
for dest in dist_df.keys():
if src == dest: continue
if optimum_var_dict[src, dest] != 0:
optimal_dict[(src, dest)] = 1
display(model.getAttr('x', cont_var_dict), optimal_dict)
{'Class': 5.0,
'Gym': 2.0,
'Library': 4.0,
'Pool': 2.9999999999999942,
'Ramen': 1.0}
{('Class', 'Office'): 1,
('Gym', 'Pool'): 1,
('Library', 'Class'): 1,
('Office', 'Ramen'): 1,
('Pool', 'Library'): 1,
('Ramen', 'Gym'): 1}
可視化¶
In [7]:
fig = gmaps.figure(center=(35.706219, 139.712076), zoom_level=16)
marker_drawing = gmaps.symbol_layer(locations=latlng_df, fill_color='green', stroke_color='green',
scale=5, hover_text=latlng_df.index.values,
info_box_content=latlng_df.index.values)
fig.add_layer(marker_drawing)
lines = []
for key, src_latlng in latlng_df.iterrows():
for key2, dest_latlng in latlng_df.iterrows():
if src_latlng['lat'] == dest_latlng['lat'] and src_latlng['lng'] == dest_latlng['lng']: continue
line = gmaps.Line(start=(src_latlng['lat'], src_latlng['lng']),
end=(dest_latlng['lat'], dest_latlng['lng']))
if line not in lines:
lines.append(line)
line_drawing = gmaps.drawing_layer()
line_drawing.features = lines
fig.add_layer(line_drawing)
#fig
embed_minimal_html('../_static/html/tsp_candidate.html', views=[fig])
In [8]:
%%HTML
<iframe src='../_static/html/tsp_candidate.html' width=700, height=500></iframe>
In [9]:
fig = gmaps.figure(center=(35.706219, 139.712076), zoom_level=16)
marker_drawing = gmaps.symbol_layer(locations=latlng_df, fill_color='green', stroke_color='green',
scale=5, hover_text=latlng_df.index.values, info_box_content=latlng_df.index.values)
fig.add_layer(marker_drawing)
lines = []
for key, value in optimal_dict.items():
lines.append(gmaps.Line(start=(latlng_df.loc[key[0], 'lat'], latlng_df.loc[key[0], 'lng']),
end=(latlng_df.loc[key[1], 'lat'], latlng_df.loc[key[1], 'lng'])))
line_drawing = gmaps.drawing_layer()
line_drawing.features = lines
fig.add_layer(line_drawing)
#fig
embed_minimal_html('../_static/html/tsp_optimum.html', views=[fig])
In [10]:
%%HTML
<iframe src='../_static/html/tsp_optimum.html' width=700, height=500></iframe>