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 very 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/

I wondered if I could design a technique for combinatorial optimization that uses an evolutionary algorithm. After a lot of work, I got a demo up and running.

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 were several rather tricky details. The most difficult sub-problem was finding a way to combine two permutation solutions to create a child solution.

I created a demo for the TSP with n = 20 cities. I set up a city distances mapping where distance(A, B) = 1 if A and B are adjacent and A is less than B, and distance(B, A) = 1.5 if A is greater than B. Therefore, the optimal 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 [0 1 2 5 14 4 3 6 8 18 19 16 17 15 13 12 11 10 9 7] with total distance of 1 + 1 + 3 + 9 + 15.0 + 1.5 + 3 + 2 + 10 + 1 + 4.5 + 1 + 3 + 3 + 1.5 + 1.5 + 1.5 + 1.5 + 3.0 = 67.0.

Playing cards are all about permutations and combinatorial optimization. When I was a college student at UC Irvine, my roommate Ed Koolish and I often drove to Las Vegas to play Blackjack — we knew how to card-count. Here’s an image from a project that uses machine learning to recognize cards.

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

using System;
namespace TSP_Evolutionary
{
  internal class TSP_EvoProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin evolutionary optimization ");
      Console.WriteLine("Setting TSP n = 20 cities ");
      Console.WriteLine("20! = 2,432,902,008,176,640,000 ");
      Console.WriteLine("Optimal soln is 0 1 2 . . 19 dist = 19");
      
      int numCities = 20;
      int numPop = 6;

      double[][] distances = GetDistances(numCities);

      Solver solver = new Solver(numCities, numPop, 
        distances, seed: 99);
      Console.WriteLine("\nInitial population: ");
      solver.Show();

      Console.WriteLine("\nBegin search ");
      solver.Solve(1000);
      Console.WriteLine("Done ");

      Console.WriteLine("\nFinal population: ");
      solver.Show();

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    static double[][] GetDistances(int n)
    {
      // in non-demo, these come from a file
      double[][] result = new double[n][];
      for (int i = 0; i "lt" n; ++i)
        result[i] = new double[n];
      for (int i = 0; i "lt" n; ++i)
        for (int j = 0; j "lt" n; ++j)
          if (i "lt" j) result[i][j] = (j - i) * 1.0;
          else if (i "gt" j) result[i][j] = (i - j) * 1.5;
      return result;
    }
  } // Program

  class Solver
  {
    public Random rnd;
    public int numCities; 
    public int numPop; // should be even

    public double[][] distances;

    public int[][] pop;  // array of solns[]
    public double[] errs;

    public int[] bestSoln;
    public double bestErr;

    public Solver(int numCities, int numPop,
      double[][] distances, int seed)
    {
      this.rnd = new Random(seed);
      this.numCities = numCities;
      this.numPop = numPop;
      this.bestErr = 0.0;
      this.distances = distances;

      this.pop = new int[numPop][];  // allocate
      for (int i = 0; i "lt" numPop; ++i)
        this.pop[i] = new int[numCities];
      this.errs = new double[numPop];

      for (int i = 0; i "lt" numPop; ++i)  // init
      {
        for (int j = 0; j "lt" numCities; ++j)
        {
          this.pop[i][j] = j;  // [0, 1, 2, . . ]
        }
        this.Shuffle(this.pop[i]);  // 
        this.errs[i] = this.ComputeError(this.pop[i]);
      }

      Array.Sort(this.errs, this.pop); 

      this.bestSoln = new int[numCities];
      for (int j = 0; j "lt" this.numCities; ++j)
        this.bestSoln[j] = this.pop[0][j];
      this.bestErr = this.errs[0];
    } // ctor

    public void Show()
    {
      for (int i = 0; i "lt" this.numPop; ++i)
      {
        for (int j = 0; j "lt" this.numCities; ++j)
        {
          Console.Write(this.pop[i][j] + " ");
        }
        Console.WriteLine(" | " + 
          this.errs[i].ToString("F4"));
      }

      Console.WriteLine("-----");
      for (int j = 0; j "lt" this.numCities; ++j)
        Console.Write(this.bestSoln[j] + " ");
      Console.WriteLine(" | " + 
        this.bestErr.ToString("F4"));
    }

    public void Shuffle(int[] soln)
    {
      // Fisher-Yates algorithm
      int n = soln.Length;
      for (int i = 0; i "lt" n; ++i)
      {
        int ri = this.rnd.Next(i, n);  // random index
        int tmp = soln[ri];
        soln[ri] = soln[i];
        soln[i] = tmp;
      }
    }

    public double ComputeError(int[] soln)
    {
      double d = 0.0;  // total distance between cities
      int n = soln.Length;  // aka numCities
      for (int i = 0; i "lt" n - 1; ++i)
      {
        int fromCity = soln[i];
        int toCity = soln[i + 1];
        d += this.distances[fromCity][toCity];
      }
      return d;  
    }

    public int[] MakeChild(int idx1, int idx2)  // crossover
    {
      int[] parent1 = this.pop[idx1];
      int[] parent2 = this.pop[idx2];
      int[] result = new int[this.numCities];
      int[] used = new int[this.numCities];
      int[] candidates = new int[2 * this.numCities];

      int k = 0;
      for (int i = 0; 
        i "lt" this.numCities / 2; ++i)  // left parent1
      {
        candidates[k++] = parent1[i];
      }

      for (int i = this.numCities / 2;
        i "lt" this.numCities; ++i)  // right parent2
      {
        candidates[k++] = parent2[i];
      }

      for (int i = 0; 
        i "lt" this.numCities / 2; ++i)  // left parent2
      {
        candidates[k++] = parent2[i];
      }

      for (int i = this.numCities / 2;
        i "lt" this.numCities; ++i)  // right parent1
      {
        candidates[k++] = parent1[i];
      }

      k = 0;
      for (int i = 0; i "lt" this.numCities; ++i)
      {
        while (used[candidates[k]] == 1) 
          k += 1;
        result[i] = candidates[k];
        used[candidates[k]] = 1;
      }

      return result;
    } // MakeChild

    public void Mutate(int[] soln)
    {
      int idx = this.rnd.Next(0, this.numCities - 1);
      int tmp = soln[idx];
      soln[idx] = soln[idx+1];
      soln[idx+1] = tmp;
    }

    public void Solve(int maxGen)
    {
      for (int gen = 0; gen "lt" maxGen; ++gen)
      {
        // 1. pick parent indexes
        int idx1, idx2;
        int flip = this.rnd.Next(2);
        if (flip == 0)
        {

          idx1 = this.rnd.Next(0, this.numPop / 2);
          idx2 = this.rnd.Next(this.numPop / 2, this.numPop);
        }
        else
        {
          idx2 = this.rnd.Next(0, this.numPop / 2);
          idx1 = this.rnd.Next(this.numPop / 2, this.numPop);
        }

        // 2. create a child
        int[] child = MakeChild(idx1, idx2);

        // 3. mutate
        Mutate(child);
        double childErr = this.ComputeError(child);
        
        // 4. found new best?
        if (childErr "lt" this.bestErr)
        {
          Console.WriteLine("New best soltn found at gen " +
            gen);
          for (int i = 0; i "lt" child.Length; ++i)
            this.bestSoln[i] = child[i];
          this.bestErr = childErr;
        }

        // 4. replace weak soln
        int idx = this.rnd.Next(this.numPop / 2, this.numPop);
        this.pop[idx] = child;
        this.errs[idx] = childErr;

        // 5. create an immigrant
        int[] imm = new int[this.numCities];
        for (int i = 0; i "lt" this.numCities; ++i)
          imm[i] = i;
        this.Shuffle(imm);
        double immErr = this.ComputeError(imm);

        // found new best?
        if (immErr "lt" this.bestErr)
        {
          Console.WriteLine("New best (immigrant) at gen " +
            gen);
          for (int i = 0; i "lt" imm.Length; ++i)
            this.bestSoln[i] = imm[i];
          this.bestErr = immErr;
        }

        // 4. replace weak soln
        idx = this.rnd.Next(this.numPop / 2, this.numPop);
        this.pop[idx] = imm;
        this.errs[idx] = immErr;

        // 5. sort the new population
        Array.Sort(this.errs, this.pop);

      } // gen
    } // Solve

  } // Solver

} // ns