# 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 $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\}$.

### 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 in order for $x_{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 $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.

## $n$-city TSP.​

\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 $d_{i,j}$ is distance between city $i$ and city $j$, and $x_{i,t}$ represents the $t$-th tour of city $i$ when it is 1.

## Modeling by JijModeling​

import jijmodeling as jm# define variablesd = jm.Placeholder('d', dim=2)N = d.shape.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 problemproblem = 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
\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 gcimport numpy as np# set the name list of traveling pointspoints = ['茨城県', '栃木県', '群馬県', '埼玉県', '千葉県', '東京都', '神奈川県']# get the latitude and longitudelatlng_list = []for point in points:    location = gc.osm(point)    latlng_list.append(location.latlng)# make distance matrixnum_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 matrixinst_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 samplersampler = jz.JijSASampler()# solve problemmultipliers = {"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 solutionsmax_show_num = 5print("Energy: ", results.evaluation.energy[:max_show_num])       # Energy is objective value of QUBOprint("Objective: ", results.evaluation.objective[:max_show_num]) # Objective values of original constrained problemprint("one-city violation: ", results.evaluation.constraint_violations["one-city"][:max_show_num])  # violation of constraintsprint("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 indexonecity_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 objectiveobjective = 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 solutionnonzero_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_indicestour[t_value] = i_valuetour = np.append(tour, [tour])tour
array([1, 0, 4, 5, 6, 3, 2, 1])
import matplotlib.pyplot as pltposition = np.array(geo_data['latlng_list']).Tplt.plot(*position, "o")plt.plot(position[tour], position[tour], "-")plt.show()print(np.array(geo_data['points'])[tour]) ['栃木県' '茨城県' '千葉県' '東京都' '神奈川県' '埼玉県' '群馬県' '栃木県']