Computation: Python, SAT + NN, permutation

• Imports
  

  import numpy as np
  import scipy as sp
  from ortools.sat.python import cp_model as cp
  

• Input
  

  The distance matrix

  \begin{equation*} D \in \mathbb{Z}^{n \times n}, \end{equation*}

  constructed as follows:

  1. load the instance pr76 from TSPLIB (Reinelt, n.d.)
  2. measure the Euclidean distances between all cities
  3. round the distances and convert them to integers

  n = len(PR76)
  c = np.array(PR76)[:, 1:]
  d_float = sp.spatial.distance_matrix(c, c)
  d = np.rint(d_float).astype(int)
  

  Table 1: Instance PR76 from TSPLIB Reinelt (n.d.)

  City	x	y
  1	3600	2300
  2	3100	3300
  3	4700	5750
  4	5400	5750
  5	5608	7103
  6	4493	7102
  7	3600	6950
  8	3100	7250
  9	4700	8450
  10	5400	8450
  11	5610	10053
  12	4492	10052
  13	3600	10800
  14	3100	10950
  15	4700	11650
  16	5400	11650
  17	6650	10800
  18	7300	10950
  19	7300	7250
  20	6650	6950
  21	7300	3300
  22	6650	2300
  23	5400	1600
  24	8350	2300
  25	7850	3300
  26	9450	5750
  27	10150	5750
  28	10358	7103
  29	9243	7102
  30	8350	6950
  31	7850	7250
  32	9450	8450
  33	10150	8450
  34	10360	10053
  35	9242	10052
  36	8350	10800
  37	7850	10950
  38	9450	11650
  39	10150	11650
  40	11400	10800
  41	12050	10950
  42	12050	7250
  43	11400	6950
  44	12050	3300
  45	11400	2300
  46	10150	1600
  47	13100	2300
  48	12600	3300
  49	14200	5750
  50	14900	5750
  51	15108	7103
  52	13993	7102
  53	13100	6950
  54	12600	7250
  55	14200	8450
  56	14900	8450
  57	15110	10053
  58	13992	10052
  59	13100	10800
  60	12600	10950
  61	14200	11650
  62	14900	11650
  63	16150	10800
  64	16800	10950
  65	16800	7250
  66	16150	6950
  67	16800	3300
  68	16150	2300
  69	14900	1600
  70	19800	800
  71	19800	10000
  72	19800	11900
  73	19800	12200
  74	200	12200
  75	200	1100
  76	200	800

• Model
  

  model = cp.CpModel()
  

• Decision variables
  

  A permutation of the cities, forming a tour,

  \begin{equation*} x \in \mathbb{N}^{n} \end{equation*}

  x = [model.new_int_var(0, n - 1, f"x_{i}") for i in range(n)]
  

• Auxiliary variables: Subscripts
  

  # Create auxiliary variables for subscripts with flattened distance indexes.
  s = [model.new_int_var(0, (n * n) - 1, f"s_{i}") for i in range(n)]
  
  # Constrain the values of the subscripts.
  for i in range(n):
      j = i + 1 if i < n - 1 else 0 # close the cycle
      model.add(s[i] == x[i] * n + x[j])
  

• Auxiliary variables: Terms
  

  # Create auxiliary variables for terms.
  t = [model.new_int_var(0, d.max(), f"d_{i}") for i in range(n)]
  
  # Constrain the values of the terms.
  d_flat = d.flatten().tolist()
  for i in range(n):
      model.add_element(s[i], d_flat, t[i])
  

• Hint
  

  x_ = [22, 21, 20, 24, 23, 45, 44, 43, 47, 46, 68, 67, 66, 49, 48, 51, 52, 53, 41, 42, 27, 28, 29, 30
  , 18, 19, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 36, 35, 34, 33, 39, 40, 59, 58
  , 57, 56, 62, 63, 61, 60, 54, 55, 50, 65, 64, 70, 71, 72, 38, 37, 31, 32, 26, 25, 3, 2, 1, 0
  , 74, 75, 73, 69]
  
  for i in range(n):
      j = (i + 1) % n
      model.add_hint(x[i], x_[i])
      model.add_hint(s[i], x_[i] * n + x_[j])
      model.add_hint(t[i], d[x_[i]][x_[j]])
  

  None
  

• Constraints: Permutation
  

  Each city must be visited exactly once.

  model.add_all_different(x)
  

• Constraints: Rotational symmetry
  

  Fix the first city to eliminate all rotations.

  model.add(x[0] == x_[0])
  

  N.B. There is also a reversal symmetry that is harder to define.

• Objective
  

  Minimize the tour length.

  \begin{align*} y = \sum_{i = 1}^{n - 1} d_{x_i, x_{i + 1}} + d_{x_n, x_1} \end{align*}

  N.B.

  • Indexing d with x[i] must be done with add_element.
  • Indexing d with add_element requires flattening of d.

  # Set the objective function.
  model.minimize(sum(t))
  

• Solution
  

  Find a solution in the given time.

  solver = cp.CpSolver()
  solver.parameters.max_time_in_seconds = 60 * 1
  result = solver.solve(model)
  y = solver.objective_value
  tour = [solver.value(x_i) for x_i in x]
  print(f"Result: {result}")
  print(f"Distance: {y}")
  print(f"Tour: {tour}")
  

  Result: 2
  Distance: 148816.0
  Tour: [22, 21, 20, 24, 23, 45, 44, 43, 47, 46, 68, 67, 66, 49, 48, 51, 52, 53, 41, 42, 27, 28, 29, 30, 18, 19, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 36, 35, 34, 33, 39, 40, 59, 58, 57, 56, 62, 63, 61, 60, 54, 55, 50, 65, 64, 70, 71, 72, 38, 37, 31, 32, 26, 25, 3, 2, 1, 0, 74, 75, 73, 69]