Traveling salesman problem in Qiskit: cannot find a solution on simulator

Question

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:

1. $0 \rightarrow 1 \rightarrow 2 \rightarrow 3 \rightarrow 0$
2. $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)

Answer 1

Qiskit has a tutorial documentation about TSP, you can find more detail at that site. As for the problem of your code, I suggest you use the qiskit-built function

tsp.random_tsp(3,seed=123) # 3 for three cities

to generate the route, instead of a single distance matrix you have written. Because tsp.random_tsp(3,seed=123) generates the coordinates and distance matrix correspondingly, while your coordinates are simply a $i$ iteration.

To see the difference between the two methods,

you can use the coordinates to calculate your density matrix, while your own method will be problematic at this phase.

(A list of coordinates [[1,0],[2,0],[3,0],[4,0]] cannot generate a distance matrix [[0,207,92,131],[207,0,300,350],[92,300,0,82],[131,350,82,0]] unless your tsp space is highly inhomogeneous.)

Here ends my investigation, I have not tested further. Discussions are welcomed.

Answer 2

May be this should be added as a comment, but I do not have the privilege to add comments.
Any way, I ran your code as it is and there was no errors:

Order:  (0, 1, 2, 3)  Distance:  720
Order:  (0, 1, 3, 2)  Distance:  731
Order:  (0, 2, 1, 3)  Distance:  873
Order:  (0, 2, 3, 1)  Distance:  731
Order:  (0, 3, 1, 2)  Distance:  873
Order:  (0, 3, 2, 1)  Distance:  720
Route length:  720
Route:  (0, 1, 2, 3)
Route length:  720.0
Route:  [1. 2. 3. 0.]

OS: Windows 10
Qiskit version: 0.23.0