Traveling Salesman Problem#
The Travelling Salesman Problem (TSP) is the problem of finding the shortest tour salesman to visit every city. It is known as an NP-hard problem. There are several well-known linear formulations for TSP. However, here we show a quadratic formulation for TSP because it is more suitable for Annealing.
Mathematical model#
Let’s consider \(n\)-city TSP.
Decision variable#
A binary variables \(x_{i,t}\) are defined as: $\( x_{i,t} = \begin{cases} 1~\text{the salesman visits \)t\(-th city \)i\(,}\\ 0~\text{otherwise}\\ \end{cases} \)\( for all \)i \in {0, …, n-1}\( and \)t \in {0, …, n-1}\(. \)x_{i,t}\( represents the salesman visits city \)i\( at \)t$ when it is 1.
Constraints#
We have to consider two constraints;
A city is visited only once. $\( \sum_t x_{i, t} = 1, ~\forall i \)$
The salesman visits only one city at a time. $\( \sum_i x_{i, t} = 1, ~\forall t \)$
The following figure illustrates why these two constraints are needed for \(x_{i,t}\) to represent a tour.
Objective function#
The total distance of the tour, which is an objective function to be minimized, can be written as follows using the product of \(x\) representing the edges used between \(t\) and \(t+1\).
\[
\sum_{i,j,t} d_{i, j}x_{i, t}x_{j, (t+1)\mod n}
\]
where \(d_{i, j} \geq 0\) is a distance between city \(i\) and city \(j\). “\(\mod n\)” is used to include in the sum the distance from the last city visited back to the first city visited. For example, when \(n\) is equal to 5 and node \(0\) is the first city visited, \(x_{0, (4+1)\mod 5}\) means \(x_{0, 0}\).
\(n\)-city TSP.#
\[\begin{split}
\begin{aligned}
\min_x \sum_{i, j} &d_{i,j} \sum_t x_{i,t} x_{j, (t+1) \mod n}\\
\mathrm{s.t.}~&\sum_i x_{i, t} = 1,~\forall t\\
&\sum_t x_{i, t} = 1, ~\forall i\\
&x_{i, t} \in \{0, 1\}
\end{aligned}
\end{split}\]
where \(d_{i,j}\) is distance between city \(i\) and city \(j\).
Modeling by JijModeling#
Next, we show an implementation using JijModeling. We first define variables for the mathematical model described above.
import jijmodeling as jm
# define variables
d = jm.Placeholder('d', ndim=2)
N = d.len_at(0, latex="N")
i = jm.Element('i', belong_to=(0, N))
j = jm.Element('j', belong_to=(0, N))
t = jm.Element('t', belong_to=(0, N))
x = jm.BinaryVar('x', shape=(N, N))
d is a two-dimensional array representing the distance between each city.
N denotes the number of cities. This can be obtained from the number of rows in d.
i, j and t are subscripts used in the mathematical model described above.
Finally, we define the binary variable x to be used in this optimization problem.
# set problem
problem = jm.Problem('TSP')
problem += jm.sum([i, j], d[i, j] * jm.sum(t, x[i, t]*x[j, (t+1) % N]))
problem += jm.Constraint("one-city", x[:, t].sum() == 1, forall=t)
problem += jm.Constraint("one-time", x[i, :].sum() == 1, forall=i)
jm.sum([i, j], d[i, j] * jm.sum(t, x[i, t]*x[j, (t+1) % N])) represents objective function \(\sum_{i,j,t} d_{i, j}x_{i, t}x_{j, (t+1)\mod n}\).
The following jm.Constraint("one-city", x[:, t] == 1, forall=t) represents a constraint to visit one city at one time.
x[:, t] allows us to express \(\sum_i x_{i, t}\) concisely.
We can check the implementation of the mathematical model on the Jupyter Notebook.
problem
\[\begin{split}\begin{array}{cccc}\text{Problem:} & \text{TSP} & & \\& & \min \quad \displaystyle \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} d_{i, j} \cdot \sum_{t = 0}^{N - 1} x_{i, t} \cdot x_{j, \left(t + 1\right) \bmod N} & \\\text{{s.t.}} & & & \\ & \text{one-city} & \displaystyle \sum_{\ast_{0} = 0}^{N - 1} x_{\ast_{0}, t} = 1 & \forall t \in \left\{0,\ldots,N - 1\right\} \\ & \text{one-time} & \displaystyle \sum_{\ast_{1} = 0}^{N - 1} x_{i, \ast_{1}} = 1 & \forall i \in \left\{0,\ldots,N - 1\right\} \\\text{{where}} & & & \\& x & 2\text{-dim binary variable}\\\end{array}\end{split}\]
Prepare an instance#
We prepare the number of cities and their coordinates.
Here we select a metropolitan area in Japan using geocoder.
from geopy.geocoders import Nominatim
import numpy as np
geolocator = Nominatim(user_agent='python')
# set the name list of traveling points
points = ['Ibaraki', 'Tochigi', 'Gunma', 'Saitama', 'Chiba', 'Tokyo', 'Kanagawa']
# get the latitude and longitude
latlng_list = []
for point in points:
location = geolocator.geocode(point)
if location is not None:
latlng_list.append([ location.latitude, location.longitude ])
# make distance matrix
num_points = len(points)
inst_d = np.zeros((num_points, num_points))
for i in range(num_points):
for j in range(num_points):
a = np.array(latlng_list[i])
b = np.array(latlng_list[j])
inst_d[i][j] = np.linalg.norm(a-b)
# normalize distance matrix
inst_d = (inst_d-inst_d.min()) / (inst_d.max()-inst_d.min())
geo_data = {'points': points, 'latlng_list': latlng_list}
instance_data = {'d': inst_d}
---------------------------------------------------------------------------
timeout Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
533 try:
--> 534 response = conn.getresponse()
535 except (BaseSSLError, OSError) as e:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connection.py:516, in HTTPConnection.getresponse(self)
515 # Get the response from http.client.HTTPConnection
--> 516 httplib_response = super().getresponse()
518 try:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/http/client.py:1377, in HTTPConnection.getresponse(self)
1376 try:
-> 1377 response.begin()
1378 except ConnectionError:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/http/client.py:320, in HTTPResponse.begin(self)
319 while True:
--> 320 version, status, reason = self._read_status()
321 if status != CONTINUE:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/http/client.py:281, in HTTPResponse._read_status(self)
280 def _read_status(self):
--> 281 line = str(self.fp.readline(_MAXLINE + 1), "iso-8859-1")
282 if len(line) > _MAXLINE:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/socket.py:716, in SocketIO.readinto(self, b)
715 try:
--> 716 return self._sock.recv_into(b)
717 except timeout:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/ssl.py:1275, in SSLSocket.recv_into(self, buffer, nbytes, flags)
1272 raise ValueError(
1273 "non-zero flags not allowed in calls to recv_into() on %s" %
1274 self.__class__)
-> 1275 return self.read(nbytes, buffer)
1276 else:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/ssl.py:1133, in SSLSocket.read(self, len, buffer)
1132 if buffer is not None:
-> 1133 return self._sslobj.read(len, buffer)
1134 else:
timeout: The read operation timed out
The above exception was the direct cause of the following exception:
ReadTimeoutError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
786 # Make the request on the HTTPConnection object
--> 787 response = self._make_request(
788 conn,
789 method,
790 url,
791 timeout=timeout_obj,
792 body=body,
793 headers=headers,
794 chunked=chunked,
795 retries=retries,
796 response_conn=response_conn,
797 preload_content=preload_content,
798 decode_content=decode_content,
799 **response_kw,
800 )
802 # Everything went great!
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:536, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
535 except (BaseSSLError, OSError) as e:
--> 536 self._raise_timeout(err=e, url=url, timeout_value=read_timeout)
537 raise
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:367, in HTTPConnectionPool._raise_timeout(self, err, url, timeout_value)
366 if isinstance(err, SocketTimeout):
--> 367 raise ReadTimeoutError(
368 self, url, f"Read timed out. (read timeout={timeout_value})"
369 ) from err
371 # See the above comment about EAGAIN in Python 3.
ReadTimeoutError: HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Read timed out. (read timeout=1)
The above exception was the direct cause of the following exception:
MaxRetryError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/requests/adapters.py:667, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
666 try:
--> 667 resp = conn.urlopen(
668 method=request.method,
669 url=url,
670 body=request.body,
671 headers=request.headers,
672 redirect=False,
673 assert_same_host=False,
674 preload_content=False,
675 decode_content=False,
676 retries=self.max_retries,
677 timeout=timeout,
678 chunked=chunked,
679 )
681 except (ProtocolError, OSError) as err:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:871, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
868 log.warning(
869 "Retrying (%r) after connection broken by '%r': %s", retries, err, url
870 )
--> 871 return self.urlopen(
872 method,
873 url,
874 body,
875 headers,
876 retries,
877 redirect,
878 assert_same_host,
879 timeout=timeout,
880 pool_timeout=pool_timeout,
881 release_conn=release_conn,
882 chunked=chunked,
883 body_pos=body_pos,
884 preload_content=preload_content,
885 decode_content=decode_content,
886 **response_kw,
887 )
889 # Handle redirect?
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:871, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
868 log.warning(
869 "Retrying (%r) after connection broken by '%r': %s", retries, err, url
870 )
--> 871 return self.urlopen(
872 method,
873 url,
874 body,
875 headers,
876 retries,
877 redirect,
878 assert_same_host,
879 timeout=timeout,
880 pool_timeout=pool_timeout,
881 release_conn=release_conn,
882 chunked=chunked,
883 body_pos=body_pos,
884 preload_content=preload_content,
885 decode_content=decode_content,
886 **response_kw,
887 )
889 # Handle redirect?
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/connectionpool.py:841, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
839 new_e = ProtocolError("Connection aborted.", new_e)
--> 841 retries = retries.increment(
842 method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
843 )
844 retries.sleep()
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/urllib3/util/retry.py:519, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
518 reason = error or ResponseError(cause)
--> 519 raise MaxRetryError(_pool, url, reason) from reason # type: ignore[arg-type]
521 log.debug("Incremented Retry for (url='%s'): %r", url, new_retry)
MaxRetryError: HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Max retries exceeded with url: /search?q=Tokyo&format=json&limit=1 (Caused by ReadTimeoutError("HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Read timed out. (read timeout=1)"))
During handling of the above exception, another exception occurred:
ConnectionError Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/geopy/adapters.py:482, in RequestsAdapter._request(self, url, timeout, headers)
481 try:
--> 482 resp = self.session.get(url, timeout=timeout, headers=headers)
483 except Exception as error:
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/requests/sessions.py:602, in Session.get(self, url, **kwargs)
601 kwargs.setdefault("allow_redirects", True)
--> 602 return self.request("GET", url, **kwargs)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/requests/sessions.py:589, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
588 send_kwargs.update(settings)
--> 589 resp = self.send(prep, **send_kwargs)
591 return resp
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/requests/sessions.py:703, in Session.send(self, request, **kwargs)
702 # Send the request
--> 703 r = adapter.send(request, **kwargs)
705 # Total elapsed time of the request (approximately)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/requests/adapters.py:700, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
698 raise SSLError(e, request=request)
--> 700 raise ConnectionError(e, request=request)
702 except ClosedPoolError as e:
ConnectionError: HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Max retries exceeded with url: /search?q=Tokyo&format=json&limit=1 (Caused by ReadTimeoutError("HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Read timed out. (read timeout=1)"))
During handling of the above exception, another exception occurred:
GeocoderUnavailable Traceback (most recent call last)
Cell In[4], line 12
10 latlng_list = []
11 for point in points:
---> 12 location = geolocator.geocode(point)
13 if location is not None:
14 latlng_list.append([ location.latitude, location.longitude ])
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/geopy/geocoders/nominatim.py:297, in Nominatim.geocode(self, query, exactly_one, timeout, limit, addressdetails, language, geometry, extratags, country_codes, viewbox, bounded, featuretype, namedetails)
295 logger.debug("%s.geocode: %s", self.__class__.__name__, url)
296 callback = partial(self._parse_json, exactly_one=exactly_one)
--> 297 return self._call_geocoder(url, callback, timeout=timeout)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/geopy/geocoders/base.py:368, in Geocoder._call_geocoder(self, url, callback, timeout, is_json, headers)
366 try:
367 if is_json:
--> 368 result = self.adapter.get_json(url, timeout=timeout, headers=req_headers)
369 else:
370 result = self.adapter.get_text(url, timeout=timeout, headers=req_headers)
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/geopy/adapters.py:472, in RequestsAdapter.get_json(self, url, timeout, headers)
471 def get_json(self, url, *, timeout, headers):
--> 472 resp = self._request(url, timeout=timeout, headers=headers)
473 try:
474 return resp.json()
File /opt/hostedtoolcache/Python/3.9.22/x64/lib/python3.9/site-packages/geopy/adapters.py:494, in RequestsAdapter._request(self, url, timeout, headers)
492 raise GeocoderServiceError(message)
493 else:
--> 494 raise GeocoderUnavailable(message)
495 elif isinstance(error, requests.Timeout):
496 raise GeocoderTimedOut("Service timed out")
GeocoderUnavailable: HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Max retries exceeded with url: /search?q=Tokyo&format=json&limit=1 (Caused by ReadTimeoutError("HTTPSConnectionPool(host='nominatim.openstreetmap.org', port=443): Read timed out. (read timeout=1)"))
Solve by JijSolver#
We solve this problem using jijsolver.
import jijsolver
interpreter = jm.Interpreter(instance_data)
instance = interpreter.eval_problem(problem)
solution = jijsolver.solve(instance, time_limit_sec=1.0)
Visualize the solution#
Using the solution obtained, we visualize the order of that salesman’s visit.
import matplotlib.pyplot as plt
df = solution.decision_variables
x_indices = df[df["value"] == 1.0]["subscripts"]
tour = [index for index, _ in sorted(x_indices, key=lambda x: x[1])]
tour.append(tour[0])
position = np.array(geo_data['latlng_list']).T[[1, 0]]
plt.axes().set_aspect('equal')
plt.plot(*position, "o")
plt.plot(position[0][tour], position[1][tour], "-")
plt.show()
print(np.array(geo_data['points'])[tour])
['Gunma' 'Tochigi' 'Ibaraki' 'Chiba' 'Tokyo' 'Kanagawa' 'Saitama' 'Gunma']