I am trying to solve Traveling Salesman Problem (TSP) in Qiskit based on Qiskit Tutorial.
I used TSP for four cities described by this distance matrix:
$$ D = \begin{pmatrix} 0 & 207 & 92 & 131 \\ 207 & 0 & 300 & 350 \\ 92 & 300 & 0 & 82\\ 131 & 350 & 82& 0 \\ \end{pmatrix} $$
With brute force I found two optimal solutions:
- $0 \rightarrow 1 \rightarrow 2 \rightarrow 3 \rightarrow 0$
- $0 \rightarrow 3 \rightarrow 2 \rightarrow 1 \rightarrow 0$
A total distance is 720 for both solutions.
However, when I run the problem on qasm_simulator
with TSP
algorithm in qiskit.optimization.applications.ising
library, the returned solution is $0 \rightarrow 2 \rightarrow 3 \rightarrow 1 \rightarrow 0$ with distance 873. But according to matrix $D$, the total distance should be 731.
I can understand that the quantum solver cannot reach the optimal solution but I am rather confused by miscalculated total distance for the solution which was found.
So my questions is what wrong in my code? Just note that solution for example in Qiskit Tutorial was found correctly.
My second question is how to set TSP solver to reach an optimal solution? I would expect that since I use a simulator, there is no noise and in the end I would reach the optimal solution.
EDIT: It seems that if the code is rerun, the results are different. I reached the distance 731, user Egretta Thua even the optimal 720. However, the first city in solution should be the city no. 0 which was not the case both in my or Egretta code rerun.
Here is my code:
%matplotlib inline
# Importing standard Qiskit libraries and configuring account
from qiskit import QuantumCircuit, execute, Aer, IBMQ
from qiskit.compiler import transpile, assemble
from qiskit.tools.jupyter import *
from qiskit.visualization import *
#visualization tools
import matplotlib.pyplot as plt
import matplotlib.axes as axes
#other tool
import numpy as np
import networkx as nx
from itertools import permutations
#quadratic program
from qiskit.optimization import QuadraticProgram
#TSP libraries
from qiskit.optimization.applications.ising import tsp
from qiskit.optimization.applications.ising.common import sample_most_likely
#quantum computing optimization
from qiskit.optimization.converters import IsingToQuadraticProgram
from qiskit.aqua.algorithms import VQE, QAOA, NumPyMinimumEigensolver
from qiskit.optimization.algorithms import MinimumEigenOptimizer
#function for solving the TSP with brute force, i.e. generate all permutations and calc distances
def brute_force_tsp(w):
N = len(w)
#generate tuples with all permutation of numbers 1,2...N-1
#first index is zero but we want to start our travel in the first city (i.e. with index 0)
a = list(permutations(range(1,N)))
best_dist = 1e10 #distance at begining
for i in a: #for all permutations
distance = 0
pre_j = 0 #starting in city 0
for j in i: #for each element of a permutation
distance = distance + w[pre_j,j] #going from one city to another
pre_j = j #save previous city
distance = distance + w[pre_j,0] #going back to city 0
order = (0,) + i #route description (i is permutation, 0 at the begining - the first city)
print('Order: ', order, ' Distance: ', distance) #show solutions
if distance < best_dist:
best_dist = distance
best_order = order
print('Route length: ', best_dist)
print('Route: ', best_order)
return best_dist, best_order
#showing resulting route in graph
def show_tsp_graph(route):
n = len(route)
#showing the route in graph
G = nx.Graph() #graph
G.add_nodes_from(range(0,n)) #add nodes
#adding edges based on solution
for i in range(0,n-1):
G.add_edge(route[i], route[i+1])
G.add_edge(route[n-1], 0)
nx.draw_networkx(G) #show graph
#decoding binary output of QAOA to actual solution
def decodeQAOAresults(res):
n = int(len(res)**0.5)
results = np.zeros(n)
k = 0
for i in range(0,n): #each n elements refers to one time point i
for j in range(0,n): #in each time points there are all cities
#when x = 1 then the city j is visited in ith time point
if res[k] == 1: results[i] = j
k = k + 1
return results
def tspQuantumSolver(distances, backendName):
citiesNumber = len(distances)
coordinates = np.zeros([citiesNumber, 2])
for i in range(0, citiesNumber): coordinates[i][0] = i + 1
tspTask = tsp.TspData(name = 'TSP', dim = citiesNumber, w = distances, coord = coordinates)
isingHamiltonian, offset = tsp.get_operator(tspTask)
tspQubo = QuadraticProgram()
tspQubo.from_ising(isingHamiltonian, offset)
quantumProcessor = Aer.backends(name = backendName)[0]
qaoa = MinimumEigenOptimizer(QAOA(quantum_instance = quantumProcessor))
results = qaoa.solve(tspQubo)
print('Route length: ', results.fval)
route = decodeQAOAresults(results.x)
print('Route: ', route)
return results.fval, route
distMatrix = np.array([[0,207,92,131],
[207,0,300,350],
[92,300,0,82],
[131,350,82,0]
])
#brute force solution
lengthBrute, routeBrute = brute_force_tsp(distMatrix)
show_tsp_graph(routeBrute)
#quantum solution
lengthQuantum, routeQuantum = tspQuantumSolver(distMatrix, 'qasm_simulator')
show_tsp_graph(routeQuantum)