Skip to main content

View on GitHub

Open this notebook in GitHub to run it yourself
The “Travelling Salesman Problem” [1] refers to finding the shortest route between cities, given their relative distances. In a more general sense, given a weighted directed graph, find the shortest route along the graph that goes through all the cities, where the weights correspond to the distance between cities. For example, in the graph below, the route along 0→1→2→30\rightarrow 1\rightarrow 2\rightarrow 3 yields a total distance 3, which is the shortest: 
import networkx as nx  # noqa

nonedge = 5
graph = nx.DiGraph()
graph.add_nodes_from([0, 1, 2, 3])
graph.add_edges_from([(0, 1), (1, 2), (2, 1), (2, 3)], weight=1)
graph.add_edges_from([(0, 2), (1, 3)], weight=2)
pos = nx.planar_layout(graph)
nx.draw_networkx(graph, pos=pos)

labels = nx.get_edge_attributes(graph, "weight")
nx.draw_networkx_edge_labels(graph, pos, edge_labels=labels)
distance_matrix = nx.convert_matrix.to_numpy_array(graph, nonedge=nonedge)
output As with many real world problems, this task can be cast as a combinatorial optimization problem. This demo shows how to employ the Quantum Approximate Optimization Algorithm [2] on the Classiq platform to solve the Travelling Salesperson Problem.

​
Mathematical Formulation

First, model the problem mathematically. The input is the set of distances between the cities: this is given by a matrix ww, whose (i,j)(i,j) entry refers to the distance between city ii to city jj. The output of the model is an optimized route. Any route can be captured by a binary matrix xx that states at each step (row) which city was visited (column): xij={1at time step i the route is in city j0else\begin{aligned} x_{ij} = \begin{cases} 1 & \text{at time step } i \text{ the route is in city } j \\ 0 & \text{else} \end{cases}\\ \end{aligned} For example: x=(0100000100101000)\begin{aligned} x=\begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 \end{pmatrix} \end{aligned} means starting from city 1, going to city 3 and then to city 2, and ending at city
  1. The constrained optimization problem is defined as follows:
Find x, which minimizes the path distance - $ min⁡xi,p∈{0,1}Σi,jwi,jΣpxi,pxj,p+1\begin{aligned} \min_{x_{i, p} \in \{0, 1\}} \Sigma_{i, j} w_{i, j} \Sigma_p x_{i, p} x_{j, p + 1}\\ \end{aligned} (Note that the inner sum over $p$ is simply an indicator for whether to go from city $i$ to city $j$.) such that - each point is visited once - $ \begin{aligned} \forall i, \hspace{0.2cm} \Sigma_p x_{i, p} = 1\\ \end{aligned}
  • in each step only a single point is visited -
$ ∀p,Σixi,p=1\begin{aligned} \forall p, \hspace{0.2cm} \Sigma_i x_{i, p} = 1\\ \end{aligned} ∗∗Directedgraph:∗∗Insomecases,suchasthegraphabove,notallcitiesareconnected,anditismoresuitabletodescribetheproblemwithaweighted,directedgraph.Inthiscase,tofindtheshortestpath,assumethatunconnectedcitieshaveaninfinitedistancebetweenthem.Forexample,thegraphabovecorrespondstothismatrix: **Directed graph:** In some cases, such as the graph above, not all cities are connected, and it is more suitable to describe the problem with a weighted, directed graph. In this case, to find the shortest path, assume that unconnected cities have an infinite distance between them. For example, the graph above corresponds to this matrix: w=(∞12∞∞∞12∞1∞1∞∞∞∞)\begin{aligned} w=\begin{pmatrix} \infty & 1 & 2 & \infty \\ \infty & \infty & 1 & 2 \\ \infty & 1 & \infty & 1 \\ \infty & \infty & \infty & \infty \end{pmatrix} \end{aligned} Inpractice,choosealargeenoughweightratherthaninfinity. In practice, choose a large enough weight rather than infinity.

​
Solving with the Classiq Platform

Solve the problem with the Classiq platform using QAOA by defining a Pyomo model.

​

​
Building the Pyomo Model from a Matrix of Distances Input

import numpy as np  # noqa
import pyomo.core as pyo


## Define a function that gets the matrix of distances and returns a Pyomo model


def PyomoTSP(dis_mat: np.ndarray) -> pyo.ConcreteModel:
    model = pyo.ConcreteModel("TSP")

    assert dis_mat.shape[0] == dis_mat.shape[1], "error distance matrix is not square"

    NofCities = dis_mat.shape[0]  # total number of cities
    cities = range(NofCities)  # list of cities

    # Define the variable, which is the binary matrix x: x[i, j] = 1 indicates that point i is visited at step j
    model.x = pyo.Var(cities, cities, domain=pyo.Binary)

    # we add constraints
    @model.Constraint(cities)
    def each_step_visits_one_point_rule(model, ii):
        return sum(model.x[ii, jj] for jj in range(NofCities)) == 1

    @model.Constraint(cities)
    def each_point_visited_once_rule(model, jj):
        return sum(model.x[ii, jj] for ii in range(NofCities)) == 1

    # Define the Objective function
    def is_connected(i1: int, i2: int):
        return sum(model.x[i1, kk] * model.x[i2, kk + 1] for kk in cities[:-1])

    model.cost = pyo.Objective(
        expr=sum(
            dis_mat[i1, i2] * is_connected(i1, i2) for i1, i2 in model.x.index_set()
        )
    )

    return model

​

​
Generating a Specific Problem

Pick a specific problem: the graph introduced above: 
# Generate a graph that defines the problem
import networkx as nx

graph = nx.DiGraph()
graph.add_nodes_from([0, 1, 2, 3])
graph.add_edges_from([(0, 1), (1, 2), (2, 1), (2, 3)], weight=1)
graph.add_edges_from([(0, 2), (1, 3)], weight=2)
pos = nx.planar_layout(graph)
nx.draw_networkx(graph, pos=pos)

labels = nx.get_edge_attributes(graph, "weight")
nx.draw_networkx_edge_labels(graph, pos, edge_labels=labels);
output Convert the graph object into a matrix of distances and then generate a Pyomo model for this example: 
nonedge = 5  # this variable refers to how much we penalize for unconnected points
distance_matrix = nx.convert_matrix.to_numpy_array(graph, nonedge=nonedge)
tsp_model = PyomoTSP(distance_matrix)

​

​
Setting Up the Classiq Problem Instance

To solve the Pyomo model defined above, use the CombinatorialProblem Python class. Under the hood it translates the Pyomo model to a quantum model of QAOA [1], with the cost Hamiltonian translated from the Pyomo model. Choose the number of layers for the QAOA ansatz using the num_layers argument: 
from classiq import *
from classiq.applications.combinatorial_optimization import CombinatorialProblem

combi = CombinatorialProblem(pyo_model=tsp_model, num_layers=8)

qmod = combi.get_model()

​

​
Synthesizing the QAOA Circuit and Solving the Problem

Synthesize and view the QAOA circuit (ansatz) used to solve the optimization problem: 
qprog = combi.get_qprog()
show(qprog)
Output:

Quantum program link: https://platform.classiq.io/circuit/2zJRXpc3YfcnB9t3wHdi4C956A7
Solve the problem by calling the optimize method of the CombinatorialProblem object. For the classical optimization part of QAOA, define the maximum number of classical iterations (maxiter) and the α\alpha-parameter (quantile) for running CVaR-QAOA, an improved variation of QAOA [2]: 
optimized_params = combi.optimize(maxiter=150, quantile=0.6)
Output:

Optimization Progress:  93%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉            | 140/150 [13:14<00:56,  5.68s/it]
Check the convergence of the run: 
import matplotlib.pyplot as plt

plt.plot(combi.cost_trace)
plt.xlabel("Iterations")
plt.ylabel("Cost")
plt.title("Cost convergence")
Output:

Text(0.5, 1.0, 'Cost convergence')
output

​
Optimization Results

Examine the statistics of the algorithm. To get samples with the optimized parameters, call the sample method: 
optimization_result = combi.sample(optimized_params)
optimization_result.sort_values(by="cost").head(5)
solutionprobabilitycost
108{‘x’: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0…0.0004883
41{‘x’: [[0, 0, 0, 1], [0, 0, 0, 1], [1, 0, 0, 0…0.0009775
896{‘x’: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 0…0.0004887
147{‘x’: [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0…0.0004888
1240{‘x’: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1…0.0004888
Compare the optimized results with uniformly sampled results: 
uniform_result = combi.sample_uniform()
And compare the histograms: 
optimization_result["cost"].plot(
    kind="hist",
    bins=50,
    edgecolor="black",
    weights=optimization_result["probability"],
    alpha=0.6,
    label="optimized",
)
uniform_result["cost"].plot(
    kind="hist",
    bins=50,
    edgecolor="black",
    weights=uniform_result["probability"],
    alpha=0.6,
    label="uniform",
)
plt.legend()
plt.ylabel("Probability", fontsize=16)
plt.xlabel("cost", fontsize=16)
plt.tick_params(axis="both", labelsize=14)
output Plot the solution: 
best_solution = optimization_result.solution[optimization_result.cost.idxmin()]
best_solution
Output:
{'x': [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]}
Lastly, compare with the classical solution of the problem: 
from pyomo.opt import SolverFactory

solver = SolverFactory("couenne")
solver.solve(tsp_model)

tsp_model.display()
Output:

Model TSP

    Variables:
      x : Size=16, Index=x_index
          Key    : Lower : Value : Upper : Fixed : Stale : Domain
          (0, 0) :     0 :   1.0 :     1 : False : False : Binary
          (0, 1) :     0 :   0.0 :     1 : False : False : Binary
          (0, 2) :     0 :   0.0 :     1 : False : False : Binary
          (0, 3) :     0 :   0.0 :     1 : False : False : Binary
          (1, 0) :     0 :   0.0 :     1 : False : False : Binary
          (1, 1) :     0 :   1.0 :     1 : False : False : Binary
          (1, 2) :     0 :   0.0 :     1 : False : False : Binary
          (1, 3) :     0 :   0.0 :     1 : False : False : Binary
          (2, 0) :     0 :   0.0 :     1 : False : False : Binary
          (2, 1) :     0 :   0.0 :     1 : False : False : Binary
          (2, 2) :     0 :   1.0 :     1 : False : False : Binary
          (2, 3) :     0 :   0.0 :     1 : False : False : Binary
          (3, 0) :     0 :   0.0 :     1 : False : False : Binary
          (3, 1) :     0 :   0.0 :     1 : False : False : Binary
          (3, 2) :     0 :   0.0 :     1 : False : False : Binary
          (3, 3) :     0 :   1.0 :     1 : False : False : Binary

    Objectives:
      cost : Size=1, Index=None, Active=True
          Key  : Active : Value
          None :   True :   3.0

    Constraints:
      each_step_visits_one_point_rule : Size=4
          Key : Lower : Body : Upper
            0 :   1.0 :  1.0 :   1.0
            1 :   1.0 :  1.0 :   1.0
            2 :   1.0 :  1.0 :   1.0
            3 :   1.0 :  1.0 :   1.0
      each_point_visited_once_rule : Size=4
          Key : Lower : Body : Upper
            0 :   1.0 :  1.0 :   1.0
            1 :   1.0 :  1.0 :   1.0
            2 :   1.0 :  1.0 :   1.0
            3 :   1.0 :  1.0 :   1.0
If you get the right solution, plot it: 
best_classical_solution = np.array(
    [int(pyo.value(tsp_model.x[idx])) for idx in np.ndindex(distance_matrix.shape)]
)


best_quantum_solution = np.array(
    optimization_result.solution[optimization_result.cost.idxmin()]["x"]
)


if (best_classical_solution == best_quantum_solution.flatten()).all():
    routesol = np.zeros(4, dtype=int)
    for k in range(4):
        indices = np.where(
            best_quantum_solution.reshape(distance_matrix.shape)[k, :] == 1
        )[0]
        routesol[k] = int(indices[0]) if len(indices) > 0 else 0
    edgesol = list(nx.utils.pairwise(routesol))
    nx.draw_networkx(
        graph,
        pos,
        with_labels=True,
        edgelist=edgesol,
        edge_color="red",
        node_size=200,
        width=3,
    )
    nx.draw_networkx(graph, pos=pos)

    labels = nx.get_edge_attributes(graph, "weight")
    nx.draw_networkx_edge_labels(graph, pos, edge_labels=labels)
    print("The route of the traveller is:", routesol)
Output:

The route of the traveller is: [0 1 2 3]
output

​
References

[1] Travelling Salesman Problem (Wikipedia). [2] Farhi, Edward, Jeffrey Goldstone, and Sam Gutmann. (2014). A quantum approximate optimization algorithm. arXiv preprint arXiv:1411.4028. [3] Barkoutsos, Panagiotis Kl, et al. (2020). Improving variational quantum optimization using CVaR. Quantum 4: 256.