Rust Learning from Zero (24) —— Travelling Salesperson Problem using Evolutionary Approach with Rust

I believe that for almost anyone studying in computer science would know this famous problem below.

The travelling salesperson problem asks the following question: "Given a list of cities and the distances between each pair of cities, what is the shortest possible route that visits each city exactly once and returns to the origin city?" Wikipedia

I know that there are thousands if not millions approaches that we can use to solve this problem. But just for personal notes and demonstration, I will solve TSP using one evolutionary approach with Rust.

The idea of today's evolutionary approach (EA) is rather simple:

  1. Randomly initialise some solutions and compute cost (described as fitness in EA) respectively.
  2. Mutate one of them, compute the fitness of mutated solution. If the new fitness is worser than the last one, go step 3. Otherwise go step 4.
  3. Reject this mutation and go step 2.
  4. Accept it as a solution, discard the worst solution and repeat step 2.

For $n$ nodes in a path, there are maximum $n!$ (actually, $(n-1)!$) paths to select from. Therefore, if what we want if the global optimal, then this algorithm has to iterate all possible paths, and that gives us $O(n!)$ performance. However, EA could actually gives us a quite good suboptimal in much less time because of the way it works. One use case of EA is when we only have limited time, computation performance and we can accept suboptimal as a solution.

For example, if we have the following adjacency matrix in TSP,

ABCDE
A57415
B53410
C7327
D4429
E151079

Then we solve this problem with EA, the output could be as below.

Initial solution [0]: fitness => 33, seq => CABED
Initial solution [1]: fitness => 32, seq => CADBE
Initial solution [2]: fitness => 34, seq => CDEAB
Initial solution [3]: fitness => 32, seq => ABDEC
Initial solution [4]: fitness => 33, seq => ADEBC
[0] Reject mutation: fitness => 38, seq => DBEAC
[1] Reject mutation: fitness => 34, seq => CDAEB
[2] Accept mutation: fitness => 33, seq => CEADB
[2] Worst solution dropped: fitness => 34
[3] Accept mutation: fitness => 28, seq => EDABC
[3] Worst solution dropped: fitness => 33
[4] Reject mutation: fitness => 33, seq => CBDAE
[5] Reject mutation: fitness => 38, seq => ACDBE
Best solution so far: fitness => 28, seq => EDABC

As you can see, EA found quite a good solution with fitness 28, which is EDABC. Also, this is in fact the global optimal and EA actually found the global optimal in the 3rd mutation!

The Rust code I used is shown below.

use std::collections::HashMap;
use std::vec::Vec;
use rand::seq::SliceRandom;
use rand::thread_rng;
use std::string::String;
use std::sync::Arc;
use std::cmp::Ordering;
use factorial::Factorial;

fn main() {
    // create a 5x5 2D matrix to store path cost from one vertex to another
    let num_node = 5;
    let mut matrix: Vec<Vec<i32>> = vec![vec![0i32; num_node]; num_node];
    add_edge(&mut matrix, 0, 1, 5);
    add_edge(&mut matrix, 0, 1, 5);
    add_edge(&mut matrix, 0, 2, 7);
    add_edge(&mut matrix, 0, 3, 4);
    add_edge(&mut matrix, 0, 4, 15);
    add_edge(&mut matrix, 1, 2, 3);
    add_edge(&mut matrix, 1, 3, 4);
    add_edge(&mut matrix, 1, 4, 10);
    add_edge(&mut matrix, 2, 3, 2);
    add_edge(&mut matrix, 2, 4, 7);
    add_edge(&mut matrix, 3, 4, 9);

    // get mutation functor
    let mut next_mutation = mutator((num_node as u32).clone(), matrix);
    // create random initial solutions
    let mut solutions : Vec<(i32, Vec<i32>, String)> = vec!();
    for i in 0..num_node.clone() {
        let initial_solution = next_mutation();
        println!("Initial solution [{}]: fitness => {}, seq => {}", i, initial_solution.0, initial_solution.2);
        solutions.push(initial_solution);
    }
    sort_solution_by_fitness(&mut solutions);

    // number of mutations for update
    let mut num_mutation = 6;
    // the maximum different ways in TSP would be n!
    if num_mutation > num_node.factorial() {
        num_mutation = num_node.factorial() - num_node;
    }
    // remember worst fitness among initial solutions
    let mut worst_fitness = solutions.get(num_node - 1).unwrap().0;
    for i in 0..num_mutation {
        let update_solution = next_mutation();
        if update_solution.0 >= worst_fitness {
            println!("[{}] Reject mutation: fitness => {}, seq => {}", i, update_solution.0, update_solution.2);
        } else {
            println!("[{}] Accept mutation: fitness => {}, seq => {}", i, update_solution.0, update_solution.2);
            solutions.push(update_solution);
            sort_solution_by_fitness(&mut solutions);
            worst_fitness = drop_worst_solution(&mut solutions);
            println!("[{}] Worst solution dropped: fitness => {}", i, worst_fitness);
        }
    }

    // print best solution so far
    let best = solutions.get(0).unwrap();
    println!("Best solution so far: fitness => {}, seq => {}", best.0, best.2);
}

/// Add edge to adjacency matrix
fn add_edge(matrix: &mut Vec<Vec<i32>>, a: usize, b: usize, cost: i32) {
    matrix[a.clone()][b.clone()] = cost.clone();
    matrix[b][a] = cost;
}

/// Compute fitness of given sequence with adjacency matrix
fn compute_fitness(seq: &Vec<i32>, matrix: &Vec<Vec<i32>>) -> i32 {
    if seq.len() < 2 { i32::MAX }
    else {
        let head = seq[0] as usize;
        let mut cost = 0i32;
        for i in 0..(seq.len() - 1) {
            cost += matrix[seq[i] as usize][seq[i+1] as usize];
        }
        cost += matrix[seq[seq.len() - 1] as usize][head];
        cost
    }
}

/// Mutator
///
/// @param num_node, number of nodes
/// @param matrix, adjacency matrix with path cost
/// @return mutator function,
///     which gives next mutation result (fitness, sequence, string_representation) when invoke
fn mutator(num_node: u32, matrix: Vec<Vec<i32>>) -> Box<dyn FnMut() -> (i32, Vec<i32>, String)> {
    let mut ids: HashMap<String, bool> = HashMap::new();
    let matrix_ref= Arc::new(matrix);
    Box::new(move || {
        let mut rng = thread_rng();
        let mut nodes: Vec<u8> = vec!();
        for i in 0..num_node {
            nodes.push((i + 65) as u8);
        };
        let id = loop {
            nodes.shuffle(&mut rng);
            let id = String::from_utf8(nodes.clone()).unwrap();
            if !ids.contains_key(&id) {
                ids.insert(id.clone(), true);
                break id;
            }
        };
        let mutation = nodes.into_iter().map(|n| (n - 65) as i32).collect::<Vec<i32>>();
        let fitness = compute_fitness(&mutation, &matrix_ref);
        (fitness, mutation, id)
    })
}

/// Sort solution array inplace by fitness in ascending order
///
/// @param v, the solution array
fn sort_solution_by_fitness(v: &mut Vec<(i32, Vec<i32>, String)>) {
    v.sort_by(|a, b| {
        if a.0 > b.0 {
            Ordering::Greater
        } else if a.0 == b.0 {
            Ordering::Equal
        } else {
            Ordering::Less
        }
    });
}

/// Drop solution with worst fitness
///
/// @param v, the solution array
/// @return the fitness of dropped solution
fn drop_worst_solution(v: &mut Vec<(i32, Vec<i32>, String)>) -> i32 {
    let last_index = v.len() - 1;
    let worst_fitness = v[last_index].0;
    v.remove(last_index);
    worst_fitness
}

Leave a Reply

Your email address will not be published. Required fields are marked *

fourteen + 12 =