Algoritmo genético: Problema del viajante.

El Problema del Agente Viajero (TSP por sus siglas en inglés) o Problema del viajante, responde a la siguiente pregunta: Dada una lista de ciudades y las distancias entre cada par de ellas, ¿Cuál es la ruta más corta posible que visita cada ciudad exactamente una vez y regresa a la ciudad origen? Este es un problema NP-duro dentro en la optimización combinatoria, muy importante en la investigación de operaciones y en la ciencia de la computación.

Caso práctico

En el problema se presentan N! rutas posibles, aunque se puede simplificar ya que dada una ruta nos da igual el punto de partida y esto reduce el número de rutas a examinar en un factor N quedando (N-1)! Como no importa la dirección en que se desplace el viajante, el número de rutas a examinar se reduce nuevamente en un factor 2. Por lo tanto, hay que considerar (N-1)!/2 rutas posibles.

En la práctica, para un problema del viajante con 5 ciudades hay 12 rutas diferentes y no necesitamos un ordenador para encontrar la mejor ruta, pero apenas aumentamos el número de ciudades las posibilidades crece exponencialmente (en realidad, factorialmente):

• Para 10 ciudades hay 181.440 rutas diferentes.
• Para 30 ciudades hay más de 4·10^31 rutas posibles. Un ordenador que calcule un millón de rutas por segundo necesitaría 10^18 años para resolverlo. Dicho de otra forma, si se hubiera comenzado a calcular al comienzo de la creación del universo (hace unos 13.400 millones de años) todavía no se habría terminado.

Puede comprobarse que por cada ciudad nueva que incorporemos, el número de rutas se multiplica por el factor N y crece exponencialmente (factorialmente). Por ello el problema pertenece a la clase de problemas NP-completos.

En este tutorial aprenderemos a codificar un algoritmo genético que sea capaz de resolver este problema (también trataremos que lo resuelva en el tiempo de vida que nos queda)

Problema en cuestión

Un repartidor de mercancías tiene encargos en varias ciudades, las cuales se representan A-B-C-D-E, y debe iniciar y terminar en la cuidad A. ¿Qué ruta debe seguir para que la distancia sea mínima?

Este problema se puede resolver por fuerza bruta de la siguiente manera.
Listamos todos los caminos posibles

ABCDEA          ABCEDA         ABDCEA          ABDECA
ABEDCA          ABECDA         ACBDEA          ACBEDA
ACDEBA          ACDBEA         ACEDBA          ACEBDA
ADBCEA          ADBECA         ADCBEA          ADCEBA
ADECBA          ADEBCA         AEDBCA          AEDCBA
AEBCDA          AEBDCA         AECBDA          AECDBA

Ahora eliminamos los caminos inversos ya que la distancia que hay que recorrer de una ciudad a otra es la misma para ida y vuelta. Y calculamos la distancia del recorrido.

ABDCEA = 51
ABDECA = 42
ABCDEA = 69
ABEDCA = 59
ADCBEA = 64
ADCEBA = 46
ADEBCA = 55
ADECBA = 47
ADBCEA = 47
ADBECA = 37
ACDBEA = 59
ACBDEA = 60

Como vemos el camino más óptimo es ADBECA que tiene una distancia de 37.

Codificación

Para la programación del algoritmo genético voy a utilizar un programa que se llama Octave es muy parecido al Matlab con la diferencia de que es gratis. (Los códigos son compatible entre ambos). Primero vamos a crear un archivo llamado poblacionInicial donde pondremos el siguiente código:

function [poblacion] = poblacionInicial(tamano_poblacion)

  %Creación de la población inicial solo se creara el camino entre las ciudades BCDE

  %Porque el inicio y el final es la misma cuidad A

  %

  %La ciudades se representaran por números A = 1, B = 2, C = 3, D = 4, E = 5

  %

  poblacion = zeros(tamano_poblacion, 4); %Creación del arreglo población inicializando todos sus elementos en 0

 

  for fila = 1:tamano_poblacion %Recorro las filas (individuos de la población)

    for columna = 1:4 %Recorro las columnas (ciudades)

      parar = false;

      while (parar != true) %Repetir hasta que se encuentre una ciudad que no esté en el camino que se está creando

        a = floor(4*rand+2); %Genero las ciudades aleatoriamente (numero aleatorio entre 2 y 5)

        if columns(find(poblacion(fila,:) == a)) == 0 %Verificar si la cuidad no está en el camino (No se pueden repetir las ciudades en un mismo camino)

          poblacion(fila, columna) = a; %Si no está en el camino introducir la cuidad

          parar = true; %Para salir del ciclo while y generar la otra cuidad 

        end %endif

      end %endwhile

    end %endfor columna

  end %endfor fila

end %endfunction

Con este código ya hemos codificado el primer proceso Generar población inicial. Ahora continuamos con el segundo proceso Evaluar individuos que es la evaluación de los individuos, para eso vamos a crear otro fichero llamado funcionObjetivo donde pondremos el siguiente código:

function [funcion_objetivo] = funcionObjetivo (poblacion, distancias, tamano_poblacion)

  %Evaluación de todos los caminos generados. Aquí calcularemos la distancia de cada camino

 

  funcion_objetivo = zeros(tamano_poblacion, 1); %Creación del arreglo que tendrá la distancia para cada camino (individuo)

 

  p = [ones(tamano_poblacion, 1) poblacion ones(tamano_poblacion, 1)]; %Se le agregan a los camino la ciudad A(1) al principio y al final

 

  for fila = 1:tamano_poblacion %Se recorren todas las filas (caminos)

    for columna = 2:6 %Se recorren a partir de la 2 columna (ciudades)

      funcion_objetivo(fila, 1) += distancias(p(fila, columna-1), p(fila, columna)); %Se van acumulando las distancia entre ciudades

    end %endfor columna

  end %endfor fila

end %endfunction

Para el tercer proceso Selección he optado por el método de Selección por torneo. Para ello deben crear un fichero con el nombre seleccionTorneo donde pondrán el código siguiente:

function [numero_descendientes] = seleccionTorneo (funcion_objetivo, tamano_poblacion)

 

  %Selección por torneo consiste básicamente en seleccionar k individuos

  %que participaran en un torneo el cual ganara el que tenga la mejor función de adaptabilidad

 

  funcion_adaptacion = funcion_objetivo; %Creación de la función de adaptabilidad (en este caso va a ser el valor de la función objetivo)

  k=2; %Cantidad de individuos que participaran en el torneo

  v=zeros(k,2);

  numero_descendientes = zeros(tamano_poblacion,1); %Creación del arreglo que tendrá la cantidad de descendientes de cada individuo

 

  while sum(numero_descendientes)<tamano_poblacion %Mientras no se asignen todos los descendientes

    i=1;

    while i<=k %Repetir ciclo para seleccionar los padres que entran al torneo

      a=floor((tamano_poblacion)*rand+1); %Se escoge un padre aleatoriamente

      if a~= v(:,2) %Para evitar que se repita al padre en el torneo

        v(i,:)=[funcion_adaptacion(a,1) a]; %Se guarda función de adaptación y posición del padre

        i=i+1; %Si no se repite el padre entonces que se pase a escoger un nuevo padre para el torneo

      end %endif

    end %endwhile

    [r j]=min(v(:,1)); %Se busca el mejor padre del torneo (menor distancia)

    numero_descendientes(v(j,2),1) = numero_descendientes(v(j,2),1) + 1; %Incrementar un descendiente al padre ganador del torneo  

  end %endwhile

end %endfunction

Después del proceso de selección ya conocemos cuales son los individuos más actos para sobrevivir y cuantos descendientes van a tener. Lo siguiente que tenemos que hacer es el cuarto proceso Recombinación para el cual he seleccionado la recombinación en un punto. Para ello creamos un fichero llamado recombinacion1punto y ponemos el siguiente codigo:

function [poblacion] = recombinacion1punto(numero_descendientes, poblacion, factor_probabilidad_recombinacion, tamano_poblacion)

 

  %Este método de recombinación de un punto está ligeramente modificado

  %para este problema en particular, ya que tenemos el inconveniente de que las

  %ciudades en un mismo recorrido no se pueden repetir, así que si hiciéramos

  %el método como es estuviéramos violando esta regla del problema.

  %Para este problema en particular se seleccionó el punto de recombinación

  %para el cual cada padre va a mantener los primeros genes hasta este punto

  %y después del punto se van a ordenar las ciudades que quedan como esta en

  %el otro padre.

 

  numero_descendientes_temporar = numero_descendientes;

  padres = zeros(tamano_poblacion, 4);

  fila = 0;

 

  %Ordenar el arreglo de numero de descendientes da mayor a menor y guardarlo en un arreglo personal

  while max(numero_descendientes_temporar) != 0

    [valor posicion] = max(numero_descendientes_temporar);

    for cantidad = 1:valor

      fila += 1;

      padres(fila, :) = poblacion(posicion, :);

    end %endfor

    numero_descendientes_temporar(posicion, 1) = 0;

  end %endwhile

 

  x=[];

  i=0;

  while size(padres, 1) ~= 0 %Mientras no se alcance todos los individuos

    padre1 = padres(1,:); %Seleccionar el primer padre del arreglo como el padre1

    a = floor((size(padres,1)-1)*rand+2);

    padre2 = padres(a,:); %Seleccionar el segundo padre aleatorio como padre2

   

    padres(a,:) = []; %Borrar padre2 del arreglo

    padres(1,:) = []; %Borrar padre1 del arreglo

   

    if factor_probabilidad_recombinacion>=rand %Condición para evaluar si se recombinan los padres

      punto_recombinacion = floor((3)*rand+1); %Seleccionar punto de recombinación aleatoriamente

     

      %Proceso para ordenar las ciudades

      p1 = padre1;                                   

      p2 = padre2;

      for columna = 1:punto_recombinacion

        pos1 = find(p2(1, :) == padre1(1, columna));

        if pos1 !=0

          p2(pos1) = [];

        end %endif

     

        pos2 = find(p1(1, :) == padre2(1, columna));

        if pos2 != 0

          p1(pos2) = [];

        end %endif

      end %endfor

      %Guardando los hijos de la recombinación

      i=i+1;

      x(i, :) = [padre1(1,1:punto_recombinacion) p2];

      i=i+1;

      x(i, :)= [padre2(1,1:punto_recombinacion) p1];

    else %Si los padres no se recombinan crear dos hijos con la misma información de los padres

      i=i+1;

      x(i, :) = padre1; %Crear hijo uno con la misma información que el padre1

      i=i+1;

      x(i, :)= padre2; %Crear hijo dos con la misma información que el padre2

    endif

  end %endwhile

  poblacion=x; %Los hijos creados van a ser la nueva población (Nueva generación)

end %endfunction

Solo nos falta el quinto proceso Mutación que es el encargado de cambiar la información de algún hijo aleatoriamente con el objetivo de alcanzar zonas de soluciones que la recombinación haya descartado. Este proceso tiene un factor muy bajo por lo que casi nunca se va a realizar. Para ello crearemos un fichero llamado mutación en el cual ponemos el siguiente código:

function  [poblacion]=mutacion(poblacion, factor_probabilidad_mutacion, tamano_poblacion)

  if rand <= factor_probabilidad_mutacion %Condición para saber si se va a realizar la mutación   

    descendiente=floor((tamano_poblacion)*rand+1); %Seleccionar un individuo aleatoriamente

    gen1 = floor(4*rand+1); %Seleccionar gen1 aleatoriamente

    gen2 = gen1;

    while gen2 == gen1 %Mientras se seleccione el mismo gen volver a intentar 

      gen2 = floor(4*rand+1); %Seleccionar gen2 aleatoriamente

    end %endwhile     

   

    %Pasos para intercambiar genes

    temp = poblacion(descendiente, gen1);

    poblacion(descendiente, gen1) = poblacion(descendiente, gen2);

    poblacion(descendiente, gen2) = temp;

  end %endif

end %endfunction

Después de tener codificado todos los procesos del algoritmo genético tenemos que unirlos todos en un orden específico para que el algoritmo funcione correctamente. Para ello vamos a crear un fichero llamado main donde pondremos el siguiente código:

function [camino, distancia] = main()

  clear;

  clc;

  %Matriz que contiene las distancias entre ciudades

  distancias = [0 7 9 8 20; 7 0 10 4 11; 9 10 0 15 5; 15 4 15 0 17; 20 11 5 17 0];

  %Tamaño de la población

  tamano_poblacion = 10;

  %Factor de probabilidad para la recombinación

  factor_probabilidad_recombinacion=0.9;

  %Factor de probabilidad para la mutación 

  factor_probabilidad_mutacion = 0.005;

  %Variable que va a tener el menor valor de todas las generaciones (Solución)

  incumbente = 100;

  iter = 0;

  parar = false;

  r = 0;

  %Creación del arreglo población con la población inicial elegida aleatoriamente

  [poblacion] = poblacionInicial(tamano_poblacion);

  %Creación del arreglo que contiene los valores de la función objetivo

  [funcion_objetivo] = funcionObjetivo (poblacion, distancias, tamano_poblacion);

  %Ciclo que correrá hasta que se cumplan alguna de las dos condiciones

  %Condición 1: Que el valor de la incumbente no cambie por 10 generaciones.

  %Condición 2: Que se llegue a 100 generaciones

  while parar == false

    %Proceso de selección por torneo que nos devuelve el número de descendientes de cada individuo

    [numero_descendientes] = seleccionTorneo (funcion_objetivo, tamano_poblacion);

    %Proceso de recombinación que nos devuelve la nueva población (nueva generación)

    [poblacion] = recombinacion1punto(numero_descendientes, poblacion, factor_probabilidad_recombinacion, tamano_poblacion);

    %Obtener las funciones objetivos de la nueva población 

    [funcion_objetivo] = funcionObjetivo (poblacion, distancias, tamano_poblacion);

   

    %Proceso de mutación

    [poblacion]=mutacion(poblacion, factor_probabilidad_mutacion, tamano_poblacion);

    %Obtener las funciones objetivos de la nueva población por si existió la mutación

    [funcion_objetivo] = funcionObjetivo (poblacion, distancias, tamano_poblacion);

   

    %Proceso para seleccionar la Incumbente

    [valor posicion] = min(funcion_objetivo);

  

    if incumbente > valor

       incumbente = valor;

       x = poblacion(posicion, :);

    end %endif

   

    %Condicion de parada 1

    if valor==incumbente  

      r=r+1;  

      if r == 10

        parar = true;

      end %endif

    end %endif

   

    %Condición de parada 2

    if iter==100

      parar = true;

    end %endif

   

    iter=iter+1;

  end %endwhile

  %Devolver la distancia mínima que se obtuvo y el camino

  distancia = incumbente;

  camino = x;

end %endfunction

Ya con estos ficheros hemos terminado el algoritmo genetico ahora si corremos el fichero llamado main de la siguiente manera: [camino distancia] = main obtendremos el mejor camino y la distancia del mismo.

Ahora si corremos el algoritmo varias veces veremos que no siempre va a encontrar el mejor camino posible por ese motivo cree otro fichero que lo que hace es correr el algoritmo genético varias veces para seleccionar de entre todas las soluciones la mejor. Para ello vamos crear el fichero evaluacion_viajante y pondremos el siguiente código:

%Correr el algoritmo genético 10 veces y seleccionar la mejor solución

char = ['A' 'B' 'C' 'D' 'E'];

solucion = 0;

for i=1:10

  [camino, distancia] = main();

  distancias(i, 1) = distancia;

  caminos(i, :) = camino;

endfor

[distancia pos] = min(distancias);

%Mostrar la distancia

distancia

%Mostrar el camino

camino = [char(1) char(caminos(pos, :)) char(1)]

Descargar todos los archivos del proyecto aquí.