A few days ago I implemented solution to the Traveling Salesman Problem (TSP) using an evolutionary algorithm, implemented in the C# language. As I was writing the C# code, I suspected that the same problem implemented using Python would have a few significant differences. So I decided to refactor the TSP in C# language program to Python to verify my hunch.

See https://jamesmccaffreyblog.com/2022/11/17/traveling-salesman-problem-combinatorial-optimization-using-an-evolutionary-algorithm-with-csharp/ for the C# version.

A combinatorial optimization problem is one where the goal is to place items in a correct order. The classic example is the Traveling Salesman Problem (TSP). Suppose you have n = 20 cities that are numbered 0 to 19. There is a distance between each pair of cities. You want to find the shortest path that visits all cities. For example, one possible solution might be [19, 5, 0, 2, 4, . . 11].

Combinatorial optimization problems are difficult. There are 20! = 2,432,902,008,176,640,000 possible solutions if distance(A,B) != distance(B,A). This is a huge number and so you can’t try every permutation.

The most common technique for combinatorial optimization is to use Simulated Annealing (SA). See https://jamesmccaffreyblog.com/2021/12/01/the-traveling-salesman-problem-using-simulated-annealing-in-csharp/.

After a lot of work, I was able to design and implement a solution to TSP that uses an evolutionary algorithm. In very high level pseudo-code, an evolutionary algorithm looks like:

create a population of random solutions
loop many times
  pick two existing solutions
  make a child solution
  mutate the child
  replace a poor solution with new child
end-loop
return best solution found

There are several rather tricky details. I created a demo for the TSP with n = 20 cities. I set up a city distances mapping where distance(A, B) = (B – A) * 1.0 if A is less than B, and distance(B, A) = (A – B) * 1.5 if A is greater than B. For example, dist(3, 8) = (8-3) * 1 = 5.0 and dist(8, 3) = (8-3) * 1.5 = 7.5.

With this city distance mapping, the optimal solution/permutation is [0, 1, 2, . . 19] with total distance of 1 + 1 + . . + 1 = 19.

I used a population size of 6 and looped for 1,000 iterations. The algorithm found a very good, but not optimal, solution of [ 1 15 17 18 16 14 13 12 9 4 3 0 2 5 10 11 8 6 7 19 ] with total distance of 14 + 2 + 1 + 3 + 3 + 1.5 + 1.5 + 4.5 + 7.5 + 1.5 + 4.5 + 2 + 3 + 5 + 1 + 4.5 + 3 + 1 + 12 = 75.5.

One of the main differences between the Python version and the earlier C# version is that for Python I used a List of Tuple objects where each tuple item is a soln[] array of int and an err value of type float32. For C# I used parallel arrays. In Python, it’s easy to sort a list of Tuple objects, but in C# it’s easy to sort parallel arrays.

Most people are familiar with games like poker, craps, and Blackjack. All these games are essentially applied combinatorics. I ran across a very simple game called Bac Bo. The Player and the Banker each roll two dice. The side with the higher total wins. You can also bet on a tie result which has a variable payout depending on the sum of the tie.

Demo code. Replace “lt”, “gt” with Boolean operator symbols.

# tsp_evolutionary.py

import numpy as np

# -----------------------------------------------------------

class Solver():
  def __init__(self, num_cities, num_pop, distances, seed):
    self.num_cities = num_cities
    self.num_pop = num_pop
    self.rnd = np.random.RandomState(seed)

    self.distances = distances  # by ref

    # pop is a list of Tuple (soln[], err) to allow sorting
    self.pop = [] 
    for i in range(num_pop):
      soln = np.arange(num_cities, dtype=np.int64)
      self.rnd.shuffle(soln)
      err = self.compute_error(soln)
      self.pop.append( (soln, err) )

    self.pop = sorted(self.pop, key=lambda tup: tup[1]) 

    self.best_soln = np.copy(self.pop[0][0])  # soln
    self.best_err = self.pop[0][1]            # err
 
  def show(self):
    for i in range(self.num_pop):
      soln = self.pop[i][0]
      err = self.pop[i][1]
      print(soln, end="") 
      print(" | %0.4f " % err)

    print("-----")
    print(self.best_soln, end ="")
    print(" | %0.4f " % self.best_err)

  def compute_error(self, soln):
    n = self.num_cities
    d = 0.0
    for i in range(n-1):
      from_idx = soln[i]
      to_idx = soln[i+1]
      d += self.distances[from_idx][to_idx]
    return d

  def make_child(self, idx1, idx2):
    n = self.num_cities
    parent1 = self.pop[idx1][0]  # (soln, err)
    parent2 = self.pop[idx2][0]
    result = np.zeros(n, dtype=np.int64)
    used = np.zeros(n, dtype=np.int64)
    candidates = np.zeros(2 * n, dtype=np.int64)

    k = 0  # create candidates
    for i in range(0, n // 2):  # left parent1
      candidates[k] = parent1[i]; k += 1
    for i in range(n // 2, n):  # right parent2
      candidates[k] = parent2[i]; k += 1
    for i in range(0, n // 2):  # left parent2
      candidates[k] = parent2[i]; k += 1
    for i in range(n // 2, n):  # right parent1
      candidates[k] = parent1[i]; k += 1

    k = 0
    for i in range(n):
      while used[candidates[k]] == 1:
        k += 1
      result[i] = candidates[k]
      used[candidates[k]] = 1

    return result

  def mutate(self, soln):
    idx = self.rnd.randint(0, self.num_cities - 1)
    tmp = soln[idx]
    soln[idx] = soln[idx+1]
    soln[idx+1] = tmp

  def solve(self, max_gen):
    n = self.num_pop
    for gen in range(max_gen):
      # 1. pick parent indices
      flip = self.rnd.randint(0,2)  # 0 or 1
      if flip == 0:
        idx1 = self.rnd.randint(0, n // 2)
        idx2 = self.rnd.randint(n // 2, n)
      else:
        idx2 = self.rnd.randint(0, n // 2)
        idx1 = self.rnd.randint(n // 2, n)
      
      # 2. create child
      child = self.make_child(idx1, idx2)
      self.mutate(child)
      child_err = self.compute_error(child)

      # 3. is child new best?
      if child_err "lt" self.best_err:
        print("New best soltn found at gen " + str(gen))
        self.best_soln = np.copy(child)
        self.best_err = child_err

      # 4. replace a weak soln
      idx = self.rnd.randint(n // 2, n)
      self.pop[idx] = (child, child_err)

      # 5. create immigrant
      imm = np.arange(self.num_cities, dtype=np.int64)
      self.rnd.shuffle(imm)
      imm_err = self.compute_error(imm)
      if imm_err "lt" self.best_err:
        print("New best (immigrant) at gen " + str(gen))
        self.best_soln = np.copy(imm)
        self.best_err = imm_err
 
      # 6. sort new population
      self.pop = sorted(self.pop, key=lambda tup: tup[1]) 

# -----------------------------------------------------------

def get_distances(n_cities):
  result = np.zeros((n_cities, n_cities), dtype=np.float32)
  for i in range(n_cities):
    for j in range(n_cities):
      if i "lt" j: result[i][j] = (j - i) * 1.0
      elif j "lt" i: result[i][j] = (i - j) * 1.5
  return result

# -----------------------------------------------------------

def main():
  print("\nBegin TSP using evolutionary optimization ")
  print("Setting TSP n = 20 cities ")
  print("20! = 2,432,902,008,176,640,000 ")
  print("Optimal soln is 0 1 2 . . 19 with dist = 19.0 ")
      
  num_cities = 20
  num_pop = 6
  distances = get_distances(num_cities)
    
  solver = Solver(num_cities, num_pop, distances, seed=1)
  print("\nInitial population: ")
  solver.show()

  print("\nBegin search ")
  solver.solve(1000)
  print("Done ")

  print("\nFinal population: ")
  solver.show()

  print("\nEnd demo ")

# -----------------------------------------------------------

if __name__ == "__main__":
  main()