Merge pull request #230 from alejomonbar/TSP_LP

UF: Adding the traveling salesman problem (TSP) linear programming version

src/openqaoa-core/openqaoa/problems/__init__.py

@@ -3,7 +3,7 @@
 from .minimumvertexcover import MinimumVertexCover
 from .numberpartition import NumberPartition
 from .shortestpath import ShortestPath
-from .tsp import TSP
+from .tsp import TSP, TSP_LP
 from .portfoliooptimization import PortfolioOptimization
 from .maximalindependentset import MIS
 from .binpacking import BinPacking

src/openqaoa-core/openqaoa/problems/helper_functions.py

@@ -7,7 +7,7 @@
 from .minimumvertexcover import MinimumVertexCover
 from .numberpartition import NumberPartition
 from .shortestpath import ShortestPath
-from .tsp import TSP
+from .tsp import TSP, TSP_LP
 from .vehiclerouting import VRP
 from .maximalindependentset import MIS
 from .binpacking import BinPacking
@@ -37,6 +37,7 @@ def create_problem_from_dict(problem_instance: dict) -> Problem:
     problem_mapper = {
         "generic_qubo": QUBO,
         "tsp": TSP,
+        "tsp_lp": TSP_LP,
         "number_partition": NumberPartition,
         "maximum_cut": MaximumCut,
         "knapsack": Knapsack,

src/openqaoa-core/openqaoa/problems/tsp.py

@@ -4,6 +4,7 @@
 
 from ..utilities import check_kwargs
 from .problem import Problem
+from .vehiclerouting import VRP
 from .qubo import QUBO
 
 
@@ -350,3 +351,113 @@ def qubo(self):
         # Convert to Ising equivalent since variables are in {0, 1} rather than {-1, 1}
         ising_terms, ising_weights = QUBO.convert_qubo_to_ising(n, terms, weights)
         return QUBO(n, ising_terms, ising_weights, self.problem_instance)
+
+class TSP_LP(VRP):
+    """
+    Creates an instance of the traveling salesman problem (TSP) based on the 
+    linear programming (LP) formulation. Note that the LP formulation of the TSP
+    can be seen as a particular case of the Vehicle routing problem (VRP) with one vehicle.
+    https://en.wikipedia.org/wiki/Travelling_salesman_problem
+    
+
+    Parameters
+    ----------
+    G: networkx.Graph
+        Networkx graph of the problem
+    pos: list[list]
+        position x, y of each node
+    subtours: list[list]
+        if -1 (Default value): All the possible subtours are added to the constraints. Avoid it for large instances.
+        if there are subtours that want be avoided in the solution, e.g, a 8 cities
+        TSP with an optimal solution showing subtour between nodes 4, 5, and 8. It can be
+        avoided introducing the constraint subtours=[[4,5,8]]. Additional information
+        about subtours refer to https://de.wikipedia.org/wiki/Datei:TSP_short_cycles.png
+    penalty: float
+        Constraints penalty factor. If the method is 'unbalanced' three values are needed,
+        one for the equality constraints and two for the inequality constraints.
+    method: str
+        Two available methods for the inequality constraints ["slack", "unbalanced"]
+        For 'unblanced' see https://arxiv.org/abs/2211.13914
+    Returns
+    -------
+        An instance of the TSP problem for the linear programming formulation.
+    """
+    __name__ = "tsp_lp"
+
+    def __init__(self, G, pos, subtours, penalty, method):
+        super().__init__(G, pos, n_vehicles=1, subtours=subtours,
+                         method=method, penalty=penalty)
+    @staticmethod
+    def random_instance(**kwargs):
+        """
+        Creates a random instance of the traveling salesman problem, whose graph is
+        random.
+
+        Parameters
+        ----------
+        **kwargs:
+            Required keyword arguments are:
+            n_nodes: int
+                The number of nodes (vertices) in the graph.
+            method: str
+                method for the inequality constraints ['slack', 'unbalanced'].
+                For the unbalaced method refer to https://arxiv.org/abs/2211.13914
+                For the slack method: https://en.wikipedia.org/wiki/Slack_variable
+            subtours: list[list]
+                Manually add the subtours to be avoided
+            seed: int
+                Seed for the random problems.
+
+        Returns
+        -------
+            A random instance of the vehicle routing problem.
+        """
+        n_nodes = kwargs.get("n_nodes", 6)
+        seed = kwargs.get("seed", None)
+        method = kwargs.get("method", "slack")
+        if method == "slack":
+            penalty = kwargs.get("penalty", 4)
+        elif method == "unbalanced":
+            penalty = kwargs.get("penalty", [4, 1, 1])
+        else:
+            raise ValueError(f"The method '{method}' is not valid.")
+        subtours = kwargs.get("subtours", -1)
+        np.random.seed(seed)
+        G = nx.Graph()
+        G.add_nodes_from(range(n_nodes))
+        pos = [[0, 0]]
+        pos += [list(2 * np.random.rand(2) - 1) for _ in range(n_nodes - 1)]
+        for i in range(n_nodes - 1):
+            for j in range(i + 1, n_nodes):
+                r = np.sqrt((pos[i][0] - pos[j][0]) ** 2 + (pos[i][1] - pos[j][1]) ** 2)
+                G.add_weighted_edges_from([(i, j, r)])
+
+        return TSP_LP(G, pos, subtours=subtours, method=method, penalty=penalty)
+
+    def get_distance(self, sol):
+        """
+        Get the TSP solution distance
+
+        Parameters
+        ----------
+        sol : str, dict
+            Solution of the TSP problem.
+
+        Returns
+        -------
+        total_distance : float
+            Total distance of the current solution.
+
+        """
+        cities = self.G.number_of_nodes()
+        if isinstance(sol, str):
+            solution = {}
+            for n, var in enumerate(self.docplex_model.iter_binary_vars()):
+                solution[var.name] = int(sol[n])
+            sol = solution
+        total_distance = 0
+        for i in range(cities):
+            for j in range(i+1, cities):
+                if sol[f"x_{i}_{j}"]:
+                    total_distance += self.G.edges[i, j]["weight"]
+        return total_distance

src/openqaoa-core/tests/test_qubo.py

@@ -5,6 +5,7 @@
     NumberPartition,
     QUBO,
     TSP,
+    TSP_LP,
     Knapsack,
     ShortestPath,
     SlackFreeKnapsack,
@@ -24,6 +25,7 @@ class TestQUBO(unittest.TestCase):
     def __generate_random_problems(self):
         problems_random_instances = {
             "tsp": TSP.random_instance(n_cities=randint(2, 15)),
+            "tsp_lp":TSP_LP.random_instance(n_nodes=randint(2, 15)),
             "number_partition": NumberPartition.random_instance(
                 n_numbers=randint(2, 15)
             ),
@@ -66,6 +68,16 @@ def test_problem_instance(self):
 
         expected_keys = {
             "tsp": ["problem_type", "n_cities", "G", "A", "B"],
+            "tsp_lp": [
+                "problem_type",
+                "G",
+                "pos",
+                "n_vehicles",
+                "depot",
+                "subtours",
+                "method",
+                "penalty",
+            ],
             "number_partition": ["problem_type", "numbers", "n_numbers"],
             "maximum_cut": ["problem_type", "G"],
             "knapsack": [

src/openqaoa-core/tests/test_tsp.py

@@ -1,7 +1,7 @@
 import unittest
 import networkx as nx
 import numpy as np
-from openqaoa.problems import TSP
+from openqaoa.problems import TSP, TSP_LP
 
 
 def terms_list_equality(terms_list1, terms_list2):
@@ -239,6 +239,35 @@ def test_tsp_type_checking(self):
             TSP(G=G)
         self.assertEqual("Edge weights should be positive", str(e.exception))
 
+    # TESTING TSP PROBLEM CLASS
+class TestTSP_LP(unittest.TestCase):
+    """Tests for TSP LP class"""
+    def test_tsp_lp_terms_weights_constant(self):
+        """Testing TSP LP problem creation"""
+        tsp_qubo = TSP_LP.random_instance(n_nodes=3, seed=1234).qubo
+        expected_terms = [[0, 1], [0, 2], [1, 2], [0], [1], [2]]
+        expected_weights = [2.0, 2.0, 2.0, 7.66823080091817, 7.707925770071554, 7.704586691688892]
+        expected_constant = 18.919256737321383
+        self.assertTrue(terms_list_equality(expected_terms, tsp_qubo.terms))
+        self.assertEqual(expected_weights, tsp_qubo.weights)
+        self.assertEqual(expected_constant, tsp_qubo.constant)
+
+    def test_tsp_lp_length(self):
+        """Testing TSP LP problem creation"""
+        cities = 3
+        tsp = TSP_LP.random_instance(n_nodes=cities, seed=111)
+        solution = tsp.classical_solution()
+        distance_expected = 2.503342058155561
+        self.assertEqual(distance_expected, tsp.get_distance(solution))
+
+    def test_tsp_lp_plot(self):
+        """Testing TSP LP problem creation"""
+        from matplotlib.pyplot import Figure
+        cities = 6
+        tsp = TSP_LP.random_instance(n_nodes=cities, seed=123, subtours=[[1,3,4]])
+        solution = tsp.classical_solution()
+        fig = tsp.plot_solution(solution)  
+        self.assertTrue(isinstance(fig, Figure))
 
 if __name__ == "__main__":
     unittest.main()