Source code for eqc_models.sequence.tsp

# (C) Quantum Computing Inc., 2024.
"""
# Traveling Salesman Problem
## Class listing
`class TSPModel` - a base class
`class MTZTSPModel` - a concrete class for the MTZ formulation of a TSP
MTZ has been chosen for a demonstration of the integer capability of the
Dirac devices. A feasible solution will produce the valid sequence of 
visits by simply ordering the nodes by the $s_i$ variables. In the literature,
the ordering variables are named with $u$, but this is avoided because the
graph notation with $u$ being a tail node and $v$ being a head node is used
here. In conjunction with this notation, the variables $x_{ij}$ reference the
node index where the index of node $u$ is $i$ and the index of node $v$ is
$j$.
"""
from typing import (Dict, Tuple)
import numpy as np
from eqc_models.base import ConstrainedPolynomialModel, InequalitiesMixin
from eqc_models.base.operators import Polynomial
class TSPModel(ConstrainedPolynomialModel):
    """ 
    The TSPModel class implements the basics for building different formualtions of 
    TSP models.
    """
    def __init__(self, D : Dict[Tuple[int, int], float]):
        self.D = D
        self.nodes = nodes = set()
        self.edges = edges = set()
        for (u, v) in D.keys():
            nodes.add(u)
            nodes.add(v)
            assert (u, v) not in edges, "Only a single edge from tail to head is allowed"
            edges.add((u, v))
        # set N to the number of nodes
        self.N = len(nodes)
        self.variables = None
    def distance(self, i : int, j : int) -> float:
        """
        Parameters
        ----------
        i: int 
            index of first node
        Returns 
        -------
        float
        Retrieves the distance between nodes at indexes i and j. Supports asymmetric 
        (D[i, j] <> D[j, i]) distances.
        """
        return self.D[i, j]
    
    def cost(self, solution : np.ndarray) -> float:
        """ 
        Solution cost is the sum of all D values where (i,j) is chosen 
        Parameters:
            :solution: np.ndarray - An array of of 0,1 values which describe a route
                       depending on the formulation chosen.
        Returns: float
        """
        raise NotImplementedError("Subclass must implement cost method")
[docs]
class MTZTSPModel(InequalitiesMixin, TSPModel):
    """
    Using the Miller Tucker Zemlin (1960) formulation of TSP, create a model
    instance that contains variables for node order (to eliminate subtours)
    and edge choices.
    >>> D = {(1, 2): 1, (2, 1): 1, (1, 3): 2, (3, 1): 2, (2, 3): 3, (3, 2): 3}
    >>> model = MTZTSPModel(D)
    >>> model.penalty_multiplier = 10
    >>> model.distance(1, 2)
    1
    >>> model.distance(3, 2)
    3
    >>> solution = np.array([1, 0, 0, 1, 1, 0, 1, 2, 3, 2, 2, 2, 5])
    >>> model.cost(solution)
    6
    >>> lhs, rhs = model.constraints
    >>> lhs.shape
    (10, 13)
    >>> model.senses
    ['EQ', 'EQ', 'EQ', 'EQ', 'EQ', 'EQ', 'LE', 'LE', 'LE', 'LE']
    >>> rhs.shape
    (10,)
    >>> (lhs@solution - rhs == 0).all()
    True
    >>> model.upper_bound
    array([1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3])
    >>> poly = model.polynomial
    >>> poly.evaluate(solution) + model.penalty_multiplier * model.offset
    6.0
    >>> infeasible = np.array([0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 4, 2, 5])
    >>> Pl, Pq = -2 * rhs.T@lhs, lhs.T@lhs
    >>> Pl.T@solution + solution.T@Pq@solution + model.offset
    0.0
    >>> Pl.T@infeasible + infeasible.T@Pq@infeasible + model.offset > 0
    True
    >>> poly.evaluate(infeasible) + model.penalty_multiplier * model.offset > 0
    True
    >>> infeasible = np.array([1, 1, 0, 1, 0, 0, 0, 1, 2, 2, 2, 0, 3])
    >>> poly.evaluate(infeasible) + model.penalty_multiplier * model.offset > 6
    True
    >>> val1 = poly.evaluate(infeasible) + model.penalty_multiplier * model.offset
    >>> model.penalty_multiplier *= 2
    >>> poly2 = model.polynomial
    >>> val2 = poly2.evaluate(infeasible) + model.penalty_multiplier * model.offset
    >>> val1 < val2
    True
    """
    def __init__(self, D : Dict[Tuple[int, int], float]) -> None:
        super(MTZTSPModel, self).__init__(D)
        self.variables = variables = []
        coefficients = []
        indices = []
        # choose a depot node
        self.depot = depot = list(self.nodes)[0]
        self.var_ub = var_ub = []
        for (u, v) in D.keys():
            varname = f"x_{u}_{v}"
            varidx = len(variables)
            variables.append(varname)
            var_ub.append(1)
            indices.append((0, varidx+1))
            coefficients.append(D[(u, v)])
        for u in self.nodes:
            variables.append(f"s_{u}")
            var_ub.append(self.N - 1)
        self.coefficients = coefficients
        self.indices = indices
        self.max_order = 2
        # Build the constraints: Two constraints for every node, one for the
        # edge chosen to enter and another for the edge chosen to leave.
        # One constraint for every edge not leading to the depot.
        depot_in = [uv for uv in self.edges if uv[-1] == depot]
        m = 2*self.N+len(self.edges) - len(depot_in)
        n = len(self.variables)
        lhs = np.zeros((m, n), dtype=np.int64)
        rhs = np.zeros((m,), dtype=np.int64)
        senses = ["EQ" for i in range(m)]
        N = self.N
        M = int(np.sqrt(N)) #  N // 2
        # scaling these coefficients
        for idx, node in enumerate(self.nodes):
            rhs[idx] = M
            rhs[self.N + idx] = M
            for (u, v) in self.edges:
                if v == node:
                    varname = f"x_{u}_{v}"
                    varidx = self.variables.index(varname)
                    lhs[idx, varidx] = M
                elif u == node:
                    varname = f"x_{u}_{v}"
                    varidx = self.variables.index(varname)
                    lhs[self.N+idx, varidx] = M
        # build these subtour elimination constraints
        # s_i - s_j + N x_{ij} <= N - 1
        idx = 0
        for (u, v) in self.edges:
            if v != depot:
                varname = f"x_{u}_{v}"
                varidx = self.variables.index(varname)
                senses[2*self.N + idx] = "LE"
                lhs[2*self.N + idx, varidx] = self.N - 1
                vidx = self.variables.index(f"s_{v}")
                lhs[2*self.N + idx, vidx] = -1
                uidx = self.variables.index(f"s_{u}")
                lhs[2*self.N + idx, uidx] = 1
                rhs[2*self.N + idx] = self.N
                idx += 1
        # update the constraint senses
        self.senses = senses
        self.lhs = lhs
        self.rhs = rhs
[docs]
    def cost(self, solution : np.ndarray) -> float:
        """ 
        Solution cost is the sum of all D values where (i,j) is chosen 
        Parameters:
            :solution: np.ndarray - An array of of 0,1 values which describe a route
                       by setting variables representing active edges to 1.
        Returns: float
        """
        # get the route selection variables
        cost_val = 0
        for (u, v) in self.edges:
            varname = f"x_{u}_{v}"
            varidx = self.variables.index(varname)
            cost_val += solution[varidx] * self.D[(u, v)]
        return cost_val
    @property
    def upper_bound(self) -> np.ndarray:
        """ 
        For all route variables, the domain is {0, 1}. The sequence variables
        can take on values in [0, 1, 2, ..., N-1]. 
        
        """
        slacks = len([sense for sense in self.senses if sense != "EQ"])
        return np.array(self.var_ub + [self.N for i in range(slacks)])