Skip to main content

Travelling Salesman Problem

The travelling salesman problem (TSP) is the problem of finding the shortest tour salesman visit every cities.

Here we show a quadratic model for TSP. However, there are several well-known linear mathematical models that you may be interested in investigating as well.

Mathematical model

Let's consider nn-city TSP.

Decision variable

A binary variables xi,tx_{i,t} are defined as:

xi,t={1 the salesman visits t-th city i,0 otherwisex_{i,t} = \begin{cases} 1~\text{the salesman visits $t$-th city $i$,}\\ 0~\text{otherwise}\\ \end{cases}

for all i{0,...,n1}i \in \{0, ..., n-1\} and t{0,...,n1}t \in \{0, ..., n-1\}.

Constraints

We have to consider two constraints;

A city is visited only once.

txi,t=1, i\sum_t x_{i, t} = 1, ~\forall i

The salesman visits only one city at a time.

ixi,t=1, t\sum_i x_{i, t} = 1, ~\forall t

The following figure illustrates why these two constraints are needed in order for xi,tx_{i,t} to represent a tour.

Objective

The total distance of the tour, which is objective function to be minimized, can be written as follows using the product of xx representing the edges used between tt and t+1t+1.

i,j,tdi,jxi,txj,(t+1)modn\sum_{i,j,t} d_{i, j}x_{i, t}x_{j, (t+1)\mod n}

where di,j0d_{i, j} \geq 0 is a distance between city ii and city jj. "modn\mod n" is used to include in the sum the distance from the last city visited back to the first city visited.

nn-city TSP.

minxi,jdi,jtxi,txj,(t+1)modns.t. ixi,t=1, ttxi,t=1, ixi,t{0,1}\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}

where di,jd_{i,j} is distance between city ii and city jj, and xi,tx_{i,t} represents the tt-th tour of city ii when it is 1.

Modeling by JijModeling

import jijmodeling as jm

# define variables
d = jm.Placeholder('d', dim=2)
N = d.shape[0].set_latex("N")
x = jm.Binary('x', shape=(N, N))
i = jm.Element('i', (0, N))
j = jm.Element('j', (0, N))
t = jm.Element('t', (0, N))


# 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] == 1, forall=t)
problem += jm.Constraint("one-time", x[i, :] == 1, forall=i)

problem
Problem: TSPmini=0N1j=0N1di,jt=0N1xi,txj,(t+1)modNs.t.one-city:iˉ0=0N1xiˉ0,t=1,t{0,,N1}one-time:iˉ1=0N1xi,iˉ1=1,i{0,,N1}xi0,i1{0,1}\begin{alignat*}{4}\text{Problem} & \text{: TSP} \\ \min & \quad \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) \mod N} \\ \text{s.t.} & \\ & \text{one-city} :\\ &\quad \quad \sum_{ \bar{i}_{0} = 0 }^{ N - 1 } x_{\bar{i}_{0},t} = 1, \forall t \in \left\{ 0 ,\ldots , N - 1 \right\} \\[8pt] & \text{one-time} :\\ &\quad \quad \sum_{ \bar{i}_{1} = 0 }^{ N - 1 } x_{i,\bar{i}_{1}} = 1, \forall i \in \left\{ 0 ,\ldots , N - 1 \right\} \\[8pt] & x_{i_{0},i_{1}} \in \{0, 1\}\end{alignat*}

Prepare an instance

!pip install geocoder
import geocoder as gc
import numpy as np

# set the name list of traveling points
points = ['茨城県', '栃木県', '群馬県', '埼玉県', '千葉県', '東京都', '神奈川県']
# get the latitude and longitude
latlng_list = []
for point in points:
location = gc.osm(point)
latlng_list.append(location.latlng)
# 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}

Solve by JijZept's SA

import jijzept as jz

# set sampler
sampler = jz.JijSASampler()
# solve problem
multipliers = {"one-city": 0.5, "one-time": 0.5}
results = sampler.sample_model(problem, instance_data, multipliers, num_reads=100, search=True)

Check the results

  • results.record: store the value of solutions
  • results.evaluation: store the results of evaluation of the solutions.

First, check the results of evaluation.

# Show the result of evaluation of solutions
max_show_num = 5
print("Energy: ", results.evaluation.energy[:max_show_num]) # Energy is objective value of QUBO
print("Objective: ", results.evaluation.objective[:max_show_num]) # Objective values of original constrained problem
print("one-city violation: ", results.evaluation.constraint_violations["one-city"][:max_show_num]) # violation of constraints
print("one-time violation: ", results.evaluation.constraint_violations["one-time"][:max_show_num]) # violation of constraints
Energy:  [-3.95420241355896, -4.122001647949219, -4.056142807006836, -3.8330631256103516, -4.122001647949219]
Objective: [1.0457976829496218, 2.877998287861795, 1.9438573324949755, 1.1669369647653187, 2.877998287861795]
one-city violation: [3.0, 0.0, 1.0, 2.0, 0.0]
one-time violation: [1.0, 0.0, 1.0, 2.0, 0.0]

Extract feasible solutions and an index of lowest solution

import numpy as np
# Get feasible solution index
onecity_violation = np.array(results.evaluation.constraint_violations["one-city"])
onetime_violation = np.array(results.evaluation.constraint_violations["one-time"])

feasible = [i for i, violation in enumerate(onecity_violation + onetime_violation) if violation == 0]

# Get feasible objective
objective = np.array(results.evaluation.objective)
feas_obj = {i: obj_value for i, obj_value in zip(feasible, objective[feasible])}

lowest_index = min(feas_obj, key=feas_obj.get)

print(f"Lowest solution index: {lowest_index}, Lowest objective value: {feas_obj[lowest_index]}")
Lowest solution index: 9, Lowest objective value: 2.8779982878617947

Check the solution

Finally, we get the solution from JijZept.

# check solution
nonzero_indices, nonzero_values, shape = results.record.solution["x"][lowest_index]
print("indices: ", nonzero_indices)
print("values: ", nonzero_values)
indices:  ([0, 1, 2, 3, 4, 5, 6], [1, 0, 6, 5, 2, 3, 4])
values: [1, 1, 1, 1, 1, 1, 1]

Draw tour

N_value = len(geo_data['latlng_list'])
tour = np.zeros(N_value, dtype=int)

i_value, t_value = nonzero_indices
tour[t_value] = i_value

tour = np.append(tour, [tour[0]])
tour
array([1, 0, 4, 5, 6, 3, 2, 1])
import matplotlib.pyplot as plt

position = np.array(geo_data['latlng_list']).T

plt.plot(*position, "o")
plt.plot(position[0][tour], position[1][tour], "-")
plt.show()

print(np.array(geo_data['points'])[tour])

png

['栃木県' '茨城県' '千葉県' '東京都' '神奈川県' '埼玉県' '群馬県' '栃木県']