Planteamiento del problema

Optimizacion de funciones de prueba

Los métodos de optimizacion heuristica son alternativas a los metodos basados en el gradiente. Algunos, como los algoritmos evolutivos y las coloinas de hormicas, se inspiran en paradigmas de la biologia.

El objetivo de este trabajo es la aplicacion de metodos heuristicos.

Optimizacion de funciones de prueba

Funciones de prueba

Considere las siguientes funciones de prueba:

  • Función de Rosenbrock: \(f(x_1, x_2) = 100(x_1-x_2^2)+(1-x_1)^2, 2_i \in [-2.048,2.048]\), \(i?1,2\). Alcanza su valor minimo en \(x_1 = 1\) y \(x_2 = 1\).

  • Función de Rastrigin: \(f(x1,x2)= 20 + ∑_{i=1}^2 x_i^2 − 10cos(2\pi x_i) , x_i∈[−5.12,5.12]\), \(i=1,2\). Alcanza su valor minimo en \(x_1 = 0\) y \(x_2 = 0\).

  • Función de Schwefel: \(f(x1,x2) = -∑_{i=1}^2 x_i sin(\sqrt{|xi|}), xi∈[−500,500]\), \(i=1,2\). Alcanza su valor mínimo en \(x1=420.9687\) y \(x2=420.9687\).

  • Función de Griewank: \(f(x1,x2)=∑_{i=1}^2 \frac{x_i^2}{4000} − ∏_{i=1}^2 cos(\frac{x_i}{\sqrt{i}}) +1, xi∈[−600,600]\), \(i=1,2\). Alcanza su valor mínimo en \(x1=0\) y \(x2=0\).

  • Función Goldstein-Price: \(f(x_1,x_2)=[1+(x_1+x_2+1)^2(19−14x_1+3x_1^2−14x_2+6x_1x_2+3x_2^2)]\) \(×[30+(2x_1−3x_2)^2(18−32x_1+12x_1^2+48x_2−36x_1x_2+27x_2^2)]\), \(x_i∈[−2,2]\), \(i=1,2\). Alcanza su valor mínimo en \(x_1=0\) y \(x_2=−1\).

  • Funcion de las seis jorobas de cabello (six-humpcamel back): \(f(x1,x2)=(4−2.1x_1^2+x_1^{4/3})x_1^2+x_1x_2+(−4+4x_2^2)x_2^2\), \(x_1∈[−3,3]\) y \(x_2∈[−2,2]\). Alcanza su valor mínimo en \(x_1=−0.0898\) y \(x_2=0.7126\) y también en \(x_1=0.0898\) y \(x_2=0.7126\).

Optimización

  1. Escoja dos funciones de prueba

  2. Optimice las funciones en dos y tres dimensiones usando un metodo de descenso por gradiente con condicion inicial aleatoria

  3. Optimice las funciones en dos y tres dimnensiones usando: algoritmos evolutivos, optimizacion de particulas y evolucion diferencial

  4. Represente con un gif animado o un video el proceso de optimizacion de descenso por gradiente y el proceso usando el metodo heuristico.

El problema del vendedor viajero

Este este apartado se aplican metodos heuristicos para resolver el problema del vendedor viajero.

Un vendedor debe hacer un recorrido por las siguientes cuidades de colombia en su carro (no necesariamente en este orden):

  • Palmira
  • Pasto
  • Tulúa
  • Bogota
  • Pereira
  • Armenia
  • Caldas
  • Valledupar
  • Montería
  • Soledad
  • Cartagena
  • Barranquilla
  • Medellín
  • Bucaramanga
  • Cúcuta

Utilice colonias de hormigas y algotirmos geneticos para encontrar el orden optimo. El costo de desplazamiento entre ciudades es la suma del valor de la hora del vendedor (es un parámetro que debe estudiarse), el costo de los peajes y el costo del combustible. CAada equipo debe definir que carro hace el recorrido del vendedor y de allí extraer el costo del combustible.

Adicionalmente represente con un gif animado o un video cómo se comporta la mejor solución usando un grafico del recorrido del mapa de colombia.

Resultados

Optimización heurística

1.Como primer funcion de prueba se escogio la funcion de rastrigin definida a continuacion:

Función de Rastrigin: \(f(x1,x2)= 20 + ∑_{i=1}^2 x_i^2 − 10cos(2\pi x_i) , x_i∈[−5.12,5.12]\), \(i=1,2\). Alcanza su valor minimo en \(x_1 = 0\) y \(x_2 = 0\).

# Determinamos la función de rastrigin a utilizar para graficar:
f_rastrigin2d<- function(x1,x2){
  y <- 20+x1^2-10*cos(2*pi*x1)+x2^2-10*cos(2*pi*x2)
  return(y)
}

n<-50
x1 <- seq(-5.12,5.12,length.out=n)
x2 <- seq(-5.12,5.12,length.out=n)

X <- expand.grid(x1,x2)
z <- f_rastrigin2d(X[,1],X[,2])
Z <- matrix(z,ncol=n,nrow=n)

Curvas de nivel de la función de Rastrigin

#Curvas de nivel:
contour(x1,x2,Z,
        main = "Función de Rastrigin",
        sub = "Curvas de nivel de la función")

filled.contour(x1, x2, Z,
               main = "Función de Rastrigin",
              color.palette = bl2gr.colors)

#persp(x1,x2,Z,theta = 50, d=2)
#persp3D(x1,x2,Z,theta = 50, phi = 20, col.palette = bl2gr.colors)

fig <- plot_ly(z = as.matrix(Z), type = "surface")
fig <- fig %>% add_surface()
fig

Definimos la funcion vectorizada

# Rastrigin  vectorizada: para optimizar
f_rastrigin2d_vec <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  y <- 20+x1^2-10*cos(2*pi*x1)+x2^2-10*cos(2*pi*x2)
  return(y)
}

Optimizamos la funcion con evolucion diferencial, para ello utiliamos la funcion DEoptim 

# Optimizacion con evolución diferencial:
opt_ras01 <- DEoptim(fn=f_rastrigin2d_vec, lower=c(-5.12,-5.12), upper = c(5.12,5.12))
## Iteration: 1 bestvalit: 5.731015 bestmemit:   -1.052497    1.966901
## Iteration: 2 bestvalit: 5.731015 bestmemit:   -1.052497    1.966901
## Iteration: 3 bestvalit: 5.731015 bestmemit:   -1.052497    1.966901
## Iteration: 4 bestvalit: 5.731015 bestmemit:   -1.052497    1.966901
## Iteration: 5 bestvalit: 4.303786 bestmemit:   -0.030076   -0.865226
## Iteration: 6 bestvalit: 4.303786 bestmemit:   -0.030076   -0.865226
## Iteration: 7 bestvalit: 4.303786 bestmemit:   -0.030076   -0.865226
## Iteration: 8 bestvalit: 1.252860 bestmemit:    0.987321   -0.035310
## Iteration: 9 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 10 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 11 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 12 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 13 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 14 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 15 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 16 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 17 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 18 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 19 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 20 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 21 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 22 bestvalit: 0.520615 bestmemit:   -0.030076    0.041618
## Iteration: 23 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 24 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 25 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 26 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 27 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 28 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 29 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 30 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 31 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 32 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 33 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 34 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 35 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 36 bestvalit: 0.272920 bestmemit:    0.008153    0.036261
## Iteration: 37 bestvalit: 0.180107 bestmemit:   -0.030170    0.000547
## Iteration: 38 bestvalit: 0.180107 bestmemit:   -0.030170    0.000547
## Iteration: 39 bestvalit: 0.180107 bestmemit:   -0.030170    0.000547
## Iteration: 40 bestvalit: 0.180107 bestmemit:   -0.030170    0.000547
## Iteration: 41 bestvalit: 0.026855 bestmemit:   -0.011624    0.000547
## Iteration: 42 bestvalit: 0.026855 bestmemit:   -0.011624    0.000547
## Iteration: 43 bestvalit: 0.026855 bestmemit:   -0.011624    0.000547
## Iteration: 44 bestvalit: 0.026855 bestmemit:   -0.011624    0.000547
## Iteration: 45 bestvalit: 0.026855 bestmemit:   -0.011624    0.000547
## Iteration: 46 bestvalit: 0.026855 bestmemit:   -0.011624    0.000547
## Iteration: 47 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 48 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 49 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 50 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 51 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 52 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 53 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 54 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 55 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 56 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 57 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 58 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 59 bestvalit: 0.015868 bestmemit:    0.002638   -0.008546
## Iteration: 60 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 61 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 62 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 63 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 64 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 65 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 66 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 67 bestvalit: 0.005290 bestmemit:   -0.004806    0.001888
## Iteration: 68 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 69 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 70 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 71 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 72 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 73 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 74 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 75 bestvalit: 0.001820 bestmemit:    0.002352    0.001908
## Iteration: 76 bestvalit: 0.001540 bestmemit:   -0.001802    0.002126
## Iteration: 77 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 78 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 79 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 80 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 81 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 82 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 83 bestvalit: 0.000158 bestmemit:   -0.000017   -0.000893
## Iteration: 84 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 85 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 86 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 87 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 88 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 89 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 90 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 91 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 92 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 93 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 94 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 95 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 96 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 97 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 98 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 99 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 100 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 101 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 102 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 103 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 104 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 105 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 106 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 107 bestvalit: 0.000047 bestmemit:   -0.000307   -0.000381
## Iteration: 108 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 109 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 110 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 111 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 112 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 113 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 114 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 115 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 116 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 117 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 118 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 119 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 120 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 121 bestvalit: 0.000047 bestmemit:   -0.000303   -0.000381
## Iteration: 122 bestvalit: 0.000036 bestmemit:   -0.000189   -0.000381
## Iteration: 123 bestvalit: 0.000036 bestmemit:   -0.000189   -0.000381
## Iteration: 124 bestvalit: 0.000036 bestmemit:   -0.000189   -0.000381
## Iteration: 125 bestvalit: 0.000036 bestmemit:   -0.000189   -0.000381
## Iteration: 126 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 127 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 128 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 129 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 130 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 131 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 132 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 133 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 134 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 135 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 136 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 137 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 138 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 139 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 140 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 141 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 142 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 143 bestvalit: 0.000014 bestmemit:   -0.000242   -0.000112
## Iteration: 144 bestvalit: 0.000004 bestmemit:    0.000010   -0.000142
## Iteration: 145 bestvalit: 0.000004 bestmemit:    0.000010   -0.000142
## Iteration: 146 bestvalit: 0.000004 bestmemit:    0.000010   -0.000142
## Iteration: 147 bestvalit: 0.000004 bestmemit:    0.000010   -0.000142
## Iteration: 148 bestvalit: 0.000004 bestmemit:    0.000010   -0.000142
## Iteration: 149 bestvalit: 0.000004 bestmemit:    0.000010   -0.000142
## Iteration: 150 bestvalit: 0.000003 bestmemit:   -0.000086   -0.000096
## Iteration: 151 bestvalit: 0.000003 bestmemit:   -0.000086   -0.000096
## Iteration: 152 bestvalit: 0.000003 bestmemit:   -0.000116    0.000010
## Iteration: 153 bestvalit: 0.000003 bestmemit:   -0.000116    0.000010
## Iteration: 154 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 155 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 156 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 157 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 158 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 159 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 160 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 161 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 162 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 163 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 164 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 165 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 166 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 167 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 168 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 169 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 170 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 171 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 172 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 173 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 174 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 175 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 176 bestvalit: 0.000001 bestmemit:    0.000052    0.000010
## Iteration: 177 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 178 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 179 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 180 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 181 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 182 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 183 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 184 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 185 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 186 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 187 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 188 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 189 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 190 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 191 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 192 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 193 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 194 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 195 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 196 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 197 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 198 bestvalit: 0.000000 bestmemit:    0.000009    0.000010
## Iteration: 199 bestvalit: 0.000000 bestmemit:   -0.000003    0.000010
## Iteration: 200 bestvalit: 0.000000 bestmemit:   -0.000003    0.000010

Curvas de nivel de la funcion optimizada

contour(x1,x2,Z,
        main = "Función de Rastrigin",
        sub = "Curvas de nivel de la función")
lines(opt_ras01$member$pop,type="p",pch=2,col="red",lwd=3)

Optimizamos con algoritmos geneticos utilizando la funcion ga

# Modifiquemos la función de Rastrigin para minimizarla usando un método de maximización:
f_rastrigin2d_vec_inv <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  y <- 20+x1^2-10*cos(2*pi*x1)+x2^2-10*cos(2*pi*x2)
  return(-y)
}

opt_ras03 <- ga(type="real-valued", fitness = f_rastrigin2d_vec_inv, lower=c(-5.12,-5.12), upper = c(5.12,5.12), seed=5)

Resultados de la optimizacion

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  real-valued 
## Population size       =  50 
## Number of generations =  100 
## Elitism               =  2 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## Search domain = 
##          x1    x2
## lower -5.12 -5.12
## upper  5.12  5.12
## 
## GA results: 
## Iterations             = 100 
## Fitness function value = -6.90978e-06 
## Solution = 
##                x1            x2
## [1,] 0.0001823157 -3.989041e-05

Curvas de nivel

En la siguiente grafica vemos cómo el promedio, la mediana y el mejor valor de la función objetivo (el fitness) evoluciona en cada iteración repecto a la población vigente en la iteración

Optimizamos por particulas utilizando la funcion psoptim

opt_ras_pso <- psoptim(par=c(NA,NA), fn=f_rastrigin2d_vec, lower=c(-5.12,-5.12),
                             upper = c(5.12,5.12), control = list(trace.stats=TRUE,trace=1))
## S=12, K=3, p=0.2297, w0=0.7213, w1=0.7213, c.p=1.193, c.g=1.193
## v.max=NA, d=14.48, vectorize=FALSE, hybrid=off
## It 10: fitness=0.3372
## It 20: fitness=0.1666
## It 30: fitness=0.07882
## It 40: fitness=0.001879
## It 50: fitness=0.0002508
## It 60: fitness=3.737e-05
## It 70: fitness=6.386e-07
## It 80: fitness=1.302e-07
## It 90: fitness=2.267e-09
## It 100: fitness=1.17e-09
## It 110: fitness=2.172e-10
## It 120: fitness=4.245e-12
## It 130: fitness=2.894e-12
## It 140: fitness=2.565e-12
## It 150: fitness=1.776e-15
## It 160: fitness=1.776e-15
## It 170: fitness=0
## It 180: fitness=0
## It 190: fitness=0
## It 200: fitness=0
## It 210: fitness=0
## It 220: fitness=0
## It 230: fitness=0
## It 240: fitness=0
## It 250: fitness=0
## It 260: fitness=0
## It 270: fitness=0
## It 280: fitness=0
## It 290: fitness=0
## It 300: fitness=0
## It 310: fitness=0
## It 320: fitness=0
## It 330: fitness=0
## It 340: fitness=0
## It 350: fitness=0
## It 360: fitness=0
## It 370: fitness=0
## It 380: fitness=0
## It 390: fitness=0
## It 400: fitness=0
## It 410: fitness=0
## It 420: fitness=0
## It 430: fitness=0
## It 440: fitness=0
## It 450: fitness=0
## It 460: fitness=0
## It 470: fitness=0
## It 480: fitness=0
## It 490: fitness=0
## It 500: fitness=0
## It 510: fitness=0
## It 520: fitness=0
## It 530: fitness=0
## It 540: fitness=0
## It 550: fitness=0
## It 560: fitness=0
## It 570: fitness=0
## It 580: fitness=0
## It 590: fitness=0
## It 600: fitness=0
## It 610: fitness=0
## It 620: fitness=0
## It 630: fitness=0
## It 640: fitness=0
## It 650: fitness=0
## It 660: fitness=0
## It 670: fitness=0
## It 680: fitness=0
## It 690: fitness=0
## It 700: fitness=0
## It 710: fitness=0
## It 720: fitness=0
## It 730: fitness=0
## It 740: fitness=0
## It 750: fitness=0
## It 760: fitness=0
## It 770: fitness=0
## It 780: fitness=0
## It 790: fitness=0
## It 800: fitness=0
## It 810: fitness=0
## It 820: fitness=0
## It 830: fitness=0
## It 840: fitness=0
## It 850: fitness=0
## It 860: fitness=0
## It 870: fitness=0
## It 880: fitness=0
## It 890: fitness=0
## It 900: fitness=0
## It 910: fitness=0
## It 920: fitness=0
## It 930: fitness=0
## It 940: fitness=0
## It 950: fitness=0
## It 960: fitness=0
## It 970: fitness=0
## It 980: fitness=0
## It 990: fitness=0
## It 1000: fitness=0
## Maximal number of iterations reached

Curvas de nivel

Optimizamos por gradiente multivariado

# Calculamos la derivada parcial
partial_dev <- function(x,i,fun,h=0.01){
  e <- x*0 # crea un vector de ceros de la misma longitud de x
  e[i] <- h
  y <- (fun(x+e)-fun(x-e))/(2*h)
  return(y)
}

# Obtenemos cada una de las derivadas parciales de f en x:
num_grad <- function(x,fun,h=0.01){
  # x: punto del espacio donde se debe evaluar el gradiente
  # fun: función para la que se desea calcular el gradiente en x
  # h: es el tamaño de ventana para el cálculo de la derivada numérica
  d <- length(x)
  y <- mapply(FUN=partial_dev,i=1:d,MoreArgs=list(x=x,h=h,fun=fun))
  return(y)
}

Sin precondicinamiento:

# SIN PRECONDICIONAMIENTO:
#determinamos la función para evaluar el gradiente de la función Rastrigin:

deriv_grad <- function(x,fun,i=1,h=0.01){
  # x: punto en el que se evalúa el gradiente
  # fun: función para la cual se calcula la derivada del gradiente respecto a la íesima componente
  # i: i-ésima componente del vector x con respecto a la que se deriva
  e <- x*0 # crea un vector de ceros de la misma longitud de x
  e[i] <- h
  y <- (num_grad(x+e,fun=fun,h=h)-num_grad(x-e,fun=fun,h=h))/(2*h)
  return(y)
}

# Calculamos la matriz hessiana:
matriz_hessiana <- function(x,fun,h=0.01){
  # x: punto en el que se evalúa la matriz hessiana
  # fun: función a la que se le calcula la matriz hessiana en x
  # h: es el tamaño de ventana para el cálculo de la derivada numérica
  d <- length(x)
  y <- mapply(FUN=deriv_grad,i=1:d,MoreArgs=list(x=x,h=h,fun=fun),SIMPLIFY = TRUE)
  return(y)
}


# Determinamos la función de optimizador multivariado
optimizador_mult_numdev <- function(x0,fun,max_eval=100,h=0.01,eta=0.01){
  x <- matrix(NA,ncol =length(x0), nrow = max_eval)
  x[1,] <- x0
  for (i in 2:max_eval){
    num_grad_fun <- num_grad(x[i-1,],fun,h)
    H <- matriz_hessiana(x[i-1,],fun,h)
    cambio <- - eta*solve(H)%*%num_grad_fun
    x[i,] <- x[i-1,] + cambio
    cambio_opt <- sqrt(sum((x[i-1,]-x[i,])^2))
    if (cambio_opt<0.00001){
      break
    }
  }
  return(x[1:i,])
}

#pasamos rastrigin al optimizador iniciando desde (0.8,4.8)
sol_ras <- optimizador_mult_numdev(f_rastrigin2d_vec,x0=c(0.8,4.8),eta=1)

Curvas de nivel

contour(x1, x2, Z, las=1,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Función de Rastrigin",
        sub = "Curvas de nivel de la función")
lines(sol_ras, type="b",cex=1.5,col="red")

Con preacondicionamiento:

#CON PRECONDICIONAMIENTO
#hacemos el precondicionamiento para la cantidad de dimensiones
deriv_segunda <- function(x,fun,i=1,h=0.01){
  e <- x*0 # crea un vector de ceros de la misma longitud de x
  e[i] <- h
  y <- (fun(x+e)-2*fun(x)+fun(x-e))/(h^2)
  return(y)
}

# Determinamos la función de optimizador multivariado
optimizador_mult_precond <- function(x0,fun,max_eval=100,h=0.01,eta=0.01){
  x <- matrix(NA,ncol =length(x0), nrow = max_eval)
  d <- length(x0)
  x[1,] <- x0
  for (i in 2:max_eval){
    num_grad_fun <- num_grad(x[i-1,],fun,h)
    diag_H <- mapply(FUN=deriv_segunda,i=1:d,MoreArgs=list(x=x[i-1,],h=h,fun=fun))
    # print(diag_H)
    # H <- matriz_hessiana(x[i-1,],fun,h)
    # print(H)
    H_precond <- diag(1/diag_H)
    cambio <- - eta*H_precond%*%num_grad_fun
    x[i,] <- x[i-1,] + cambio
    cambio_opt <- sqrt(sum((x[i-1,]-x[i,])^2))
    if (cambio_opt<0.0000001){
      break
    }
  }
  return(x[1:i,])
}

#pasamos rastrigin al optimizador iniciando desde (0.8,4.8)
sol_ras_precon <- optimizador_mult_precond(f_rastrigin2d_vec,x0=c(0.8,4.8),eta=1)

Curvas de nivel

contour(x1, x2, Z, las=1,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Función de Rastrigin",
        sub = "Curvas de nivel de la función")
lines(sol_ras_precon, type="b",cex=1.5,col="red")
lines(sol_ras, type="b",cex=1.5,col="blue")
legend("topleft",col=c("red","blue"),legend = c("Con precondicionamiento","Sin precondicionamiento"),lty=1)

2.Como segunda funcion de prueba se escogio la funcion de Schwefel definida a continuacion:

Función de Schwefel: \(f(x1,x2) = -∑_{i=1}^2 x_i sin(\sqrt{|xi|}), xi∈[−500,500]\), \(i=1,2\). Alcanza su valor mínimo en \(x1=420.9687\) y \(x2=420.9687\).

#determinamos la función de schwefel a utilizar para graficar:
f_schwefel2d<- function(x1,x2){
  y <- (-((x1*sin(sqrt(abs(x1)))) + (x2*sin(sqrt(abs(x2)))) ))
  return(y)
}

n<-50
x1 <- seq(-500,500,length.out=n)
x2 <- seq(-500,500,length.out=n)

X <- expand.grid(x1,x2)
z <- f_schwefel2d(X[,1],X[,2])
Z <- matrix(z,ncol=n,nrow=n)

Curvas de nivel de la función de Schwefel

#Curvas de nivel:
contour(x1,x2,Z)

filled.contour(x1, x2, Z, color.palette = bl2gr.colors)

#persp(x1,x2,Z,theta = 50, d=2)
#la persp3d Esta poniendo problema para graficar con colores
#persp3d(x1,x2,Z,theta = 50, phi = 20, col.palette = bl2gr.colors)
fig <- plot_ly(z = as.matrix(Z), type = "surface")
fig <- fig %>% add_surface()
fig

Definimos la funcion vectorizada

#schwefel  vectorizada: para optimizar
f_schwefel2d_vec <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  y <- (-((x1*sin(sqrt(abs(x1)))) + (x2*sin(sqrt(abs(x2))))))
  return(y)
}

Optimizamos la funcion con evolucion diferencial, para ello utiliamos la funcion DEoptim 

# Optimizacion con evolución diferencial:
opt_schwefel_DE <- DEoptim(fn=f_schwefel2d_vec, lower=c(-500,-500), upper = c(500,500))
## Iteration: 1 bestvalit: -540.236444 bestmemit: -122.732363  423.862569
## Iteration: 2 bestvalit: -540.236444 bestmemit: -122.732363  423.862569
## Iteration: 3 bestvalit: -540.236444 bestmemit: -122.732363  423.862569
## Iteration: 4 bestvalit: -547.838472 bestmemit: -265.952057  407.877862
## Iteration: 5 bestvalit: -547.838472 bestmemit: -265.952057  407.877862
## Iteration: 6 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 7 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 8 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 9 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 10 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 11 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 12 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 13 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 14 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 15 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 16 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 17 bestvalit: -693.740293 bestmemit: -296.574524  407.877862
## Iteration: 18 bestvalit: -716.508367 bestmemit: -297.976772  419.166591
## Iteration: 19 bestvalit: -716.508367 bestmemit: -297.976772  419.166591
## Iteration: 20 bestvalit: -716.508367 bestmemit: -297.976772  419.166591
## Iteration: 21 bestvalit: -716.508367 bestmemit: -297.976772  419.166591
## Iteration: 22 bestvalit: -716.508367 bestmemit: -297.976772  419.166591
## Iteration: 23 bestvalit: -716.508367 bestmemit: -297.976772  419.166591
## Iteration: 24 bestvalit: -717.483880 bestmemit: -306.115224  419.166591
## Iteration: 25 bestvalit: -717.483880 bestmemit: -306.115224  419.166591
## Iteration: 26 bestvalit: -739.071132 bestmemit:  439.898208  441.825479
## Iteration: 27 bestvalit: -739.071132 bestmemit:  439.898208  441.825479
## Iteration: 28 bestvalit: -739.071132 bestmemit:  439.898208  441.825479
## Iteration: 29 bestvalit: -834.244725 bestmemit:  416.676750  424.301526
## Iteration: 30 bestvalit: -834.244725 bestmemit:  416.676750  424.301526
## Iteration: 31 bestvalit: -834.244725 bestmemit:  416.676750  424.301526
## Iteration: 32 bestvalit: -834.701454 bestmemit:  415.928518  421.708596
## Iteration: 33 bestvalit: -834.701454 bestmemit:  415.928518  421.708596
## Iteration: 34 bestvalit: -834.701454 bestmemit:  415.928518  421.708596
## Iteration: 35 bestvalit: -834.701454 bestmemit:  415.928518  421.708596
## Iteration: 36 bestvalit: -834.701454 bestmemit:  415.928518  421.708596
## Iteration: 37 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 38 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 39 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 40 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 41 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 42 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 43 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 44 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 45 bestvalit: -837.536335 bestmemit:  419.278171  421.708596
## Iteration: 46 bestvalit: -837.810263 bestmemit:  420.141012  421.708596
## Iteration: 47 bestvalit: -837.810263 bestmemit:  420.141012  421.708596
## Iteration: 48 bestvalit: -837.810263 bestmemit:  420.141012  421.708596
## Iteration: 49 bestvalit: -837.853831 bestmemit:  421.796825  421.417308
## Iteration: 50 bestvalit: -837.853831 bestmemit:  421.796825  421.417308
## Iteration: 51 bestvalit: -837.853831 bestmemit:  421.796825  421.417308
## Iteration: 52 bestvalit: -837.853831 bestmemit:  421.796825  421.417308
## Iteration: 53 bestvalit: -837.853831 bestmemit:  421.796825  421.417308
## Iteration: 54 bestvalit: -837.947710 bestmemit:  420.600300  421.054952
## Iteration: 55 bestvalit: -837.947710 bestmemit:  420.600300  421.054952
## Iteration: 56 bestvalit: -837.947710 bestmemit:  420.600300  421.054952
## Iteration: 57 bestvalit: -837.947710 bestmemit:  420.600300  421.054952
## Iteration: 58 bestvalit: -837.947710 bestmemit:  420.600300  421.054952
## Iteration: 59 bestvalit: -837.947710 bestmemit:  420.600300  421.054952
## Iteration: 60 bestvalit: -837.959880 bestmemit:  421.026584  420.760485
## Iteration: 61 bestvalit: -837.959880 bestmemit:  421.026584  420.760485
## Iteration: 62 bestvalit: -837.959880 bestmemit:  421.026584  420.760485
## Iteration: 63 bestvalit: -837.959880 bestmemit:  421.026584  420.760485
## Iteration: 64 bestvalit: -837.959880 bestmemit:  421.026584  420.760485
## Iteration: 65 bestvalit: -837.959880 bestmemit:  421.026584  420.760485
## Iteration: 66 bestvalit: -837.963153 bestmemit:  421.026584  420.836719
## Iteration: 67 bestvalit: -837.963153 bestmemit:  421.026584  420.836719
## Iteration: 68 bestvalit: -837.963153 bestmemit:  421.026584  420.836719
## Iteration: 69 bestvalit: -837.963153 bestmemit:  421.026584  420.836719
## Iteration: 70 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 71 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 72 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 73 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 74 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 75 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 76 bestvalit: -837.964551 bestmemit:  421.026584  420.889075
## Iteration: 77 bestvalit: -837.965388 bestmemit:  420.921667  420.997831
## Iteration: 78 bestvalit: -837.965388 bestmemit:  420.921667  420.997831
## Iteration: 79 bestvalit: -837.965388 bestmemit:  420.921667  420.997831
## Iteration: 80 bestvalit: -837.965388 bestmemit:  420.921667  420.997831
## Iteration: 81 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 82 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 83 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 84 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 85 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 86 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 87 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 88 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 89 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 90 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 91 bestvalit: -837.965760 bestmemit:  420.973150  420.978348
## Iteration: 92 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 93 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 94 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 95 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 96 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 97 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 98 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 99 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 100 bestvalit: -837.965771 bestmemit:  420.973150  420.971200
## Iteration: 101 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 102 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 103 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 104 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 105 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 106 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 107 bestvalit: -837.965772 bestmemit:  420.964107  420.968859
## Iteration: 108 bestvalit: -837.965773 bestmemit:  420.966911  420.965542
## Iteration: 109 bestvalit: -837.965773 bestmemit:  420.966911  420.965542
## Iteration: 110 bestvalit: -837.965773 bestmemit:  420.966911  420.965542
## Iteration: 111 bestvalit: -837.965773 bestmemit:  420.966911  420.965542
## Iteration: 112 bestvalit: -837.965773 bestmemit:  420.966911  420.965542
## Iteration: 113 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 114 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 115 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 116 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 117 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 118 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 119 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 120 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 121 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 122 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 123 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 124 bestvalit: -837.965774 bestmemit:  420.969328  420.969116
## Iteration: 125 bestvalit: -837.965774 bestmemit:  420.968800  420.968145
## Iteration: 126 bestvalit: -837.965774 bestmemit:  420.968800  420.968145
## Iteration: 127 bestvalit: -837.965774 bestmemit:  420.968800  420.968145
## Iteration: 128 bestvalit: -837.965774 bestmemit:  420.968800  420.968145
## Iteration: 129 bestvalit: -837.965774 bestmemit:  420.968800  420.968145
## Iteration: 130 bestvalit: -837.965775 bestmemit:  420.968800  420.968364
## Iteration: 131 bestvalit: -837.965775 bestmemit:  420.968800  420.968364
## Iteration: 132 bestvalit: -837.965775 bestmemit:  420.968921  420.968950
## Iteration: 133 bestvalit: -837.965775 bestmemit:  420.968921  420.968950
## Iteration: 134 bestvalit: -837.965775 bestmemit:  420.968704  420.968813
## Iteration: 135 bestvalit: -837.965775 bestmemit:  420.968704  420.968813
## Iteration: 136 bestvalit: -837.965775 bestmemit:  420.968722  420.968813
## Iteration: 137 bestvalit: -837.965775 bestmemit:  420.968722  420.968813
## Iteration: 138 bestvalit: -837.965775 bestmemit:  420.968722  420.968813
## Iteration: 139 bestvalit: -837.965775 bestmemit:  420.968770  420.968680
## Iteration: 140 bestvalit: -837.965775 bestmemit:  420.968713  420.968802
## Iteration: 141 bestvalit: -837.965775 bestmemit:  420.968732  420.968686
## Iteration: 142 bestvalit: -837.965775 bestmemit:  420.968722  420.968705
## Iteration: 143 bestvalit: -837.965775 bestmemit:  420.968732  420.968752
## Iteration: 144 bestvalit: -837.965775 bestmemit:  420.968732  420.968752
## Iteration: 145 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 146 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 147 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 148 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 149 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 150 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 151 bestvalit: -837.965775 bestmemit:  420.968734  420.968748
## Iteration: 152 bestvalit: -837.965775 bestmemit:  420.968750  420.968744
## Iteration: 153 bestvalit: -837.965775 bestmemit:  420.968750  420.968744
## Iteration: 154 bestvalit: -837.965775 bestmemit:  420.968750  420.968744
## Iteration: 155 bestvalit: -837.965775 bestmemit:  420.968750  420.968744
## Iteration: 156 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 157 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 158 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 159 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 160 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 161 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 162 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 163 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 164 bestvalit: -837.965775 bestmemit:  420.968745  420.968747
## Iteration: 165 bestvalit: -837.965775 bestmemit:  420.968745  420.968747
## Iteration: 166 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 167 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 168 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 169 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 170 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 171 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 172 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 173 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 174 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 175 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 176 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 177 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 178 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 179 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 180 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 181 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 182 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 183 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 184 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 185 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 186 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 187 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 188 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 189 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 190 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 191 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 192 bestvalit: -837.965775 bestmemit:  420.968746  420.968746
## Iteration: 193 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 194 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 195 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 196 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 197 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 198 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 199 bestvalit: -837.965775 bestmemit:  420.968747  420.968746
## Iteration: 200 bestvalit: -837.965775 bestmemit:  420.968747  420.968746

Curvas de nivel de la funcion optimizada

contour(x1,x2,Z,
        main = "Función de Schwefel",
        sub = "Curvas de nivel de la función")
lines(opt_schwefel_DE$member$pop,type="p",pch=2,col="red",lwd=3)

Optimizamos con algoritmos geneticos utilizando la funcion ga

#Modificamos la función de Schwefel para minimizarla usando un método de maximización:
f_schwefel2d_vec_inv <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  y <- (-((x1*sin(sqrt(abs(x1)))) + (x2*sin(sqrt(abs(x2)))) ))
  return(-y)
}

opt_schwefel_GA <- ga(type="real-valued", fitness = f_schwefel2d_vec_inv, lower=c(-500,-500), upper = c(500,500), seed=5)

Resultados de la optimizacion

## -- Genetic Algorithm ------------------- 
## 
## GA settings: 
## Type                  =  real-valued 
## Population size       =  50 
## Number of generations =  100 
## Elitism               =  2 
## Crossover probability =  0.8 
## Mutation probability  =  0.1 
## Search domain = 
##         x1   x2
## lower -500 -500
## upper  500  500
## 
## GA results: 
## Iterations             = 100 
## Fitness function value = 837.9658 
## Solution = 
##            x1       x2
## [1,] 420.9662 420.9706

Curvas de nivel

En la siguiente grafica vemos cómo el promedio, la mediana y el mejor valor de la función objetivo (el fitness) evoluciona en cada iteración repecto a la población vigente en la iteración

Optimizamos por particulas utilizando la funcion psoptim

opt_schwefel_pso <- psoptim(par=c(NA,NA), fn=f_schwefel2d_vec, lower=c(-500,-500),
                            upper = c(500,500), control = list(trace.stats=TRUE,trace=1))
## S=12, K=3, p=0.2297, w0=0.7213, w1=0.7213, c.p=1.193, c.g=1.193
## v.max=NA, d=1414, vectorize=FALSE, hybrid=off
## It 10: fitness=-717.4
## It 20: fitness=-831.2
## It 30: fitness=-832.9
## It 40: fitness=-837.9
## It 50: fitness=-837.9
## It 60: fitness=-838
## It 70: fitness=-838
## It 80: fitness=-838
## It 90: fitness=-838
## It 100: fitness=-838
## It 110: fitness=-838
## It 120: fitness=-838
## It 130: fitness=-838
## It 140: fitness=-838
## It 150: fitness=-838
## It 160: fitness=-838
## It 170: fitness=-838
## It 180: fitness=-838
## It 190: fitness=-838
## It 200: fitness=-838
## It 210: fitness=-838
## It 220: fitness=-838
## It 230: fitness=-838
## It 240: fitness=-838
## It 250: fitness=-838
## It 260: fitness=-838
## It 270: fitness=-838
## It 280: fitness=-838
## It 290: fitness=-838
## It 300: fitness=-838
## It 310: fitness=-838
## It 320: fitness=-838
## It 330: fitness=-838
## It 340: fitness=-838
## It 350: fitness=-838
## It 360: fitness=-838
## It 370: fitness=-838
## It 380: fitness=-838
## It 390: fitness=-838
## It 400: fitness=-838
## It 410: fitness=-838
## It 420: fitness=-838
## It 430: fitness=-838
## It 440: fitness=-838
## It 450: fitness=-838
## It 460: fitness=-838
## It 470: fitness=-838
## It 480: fitness=-838
## It 490: fitness=-838
## It 500: fitness=-838
## It 510: fitness=-838
## It 520: fitness=-838
## It 530: fitness=-838
## It 540: fitness=-838
## It 550: fitness=-838
## It 560: fitness=-838
## It 570: fitness=-838
## It 580: fitness=-838
## It 590: fitness=-838
## It 600: fitness=-838
## It 610: fitness=-838
## It 620: fitness=-838
## It 630: fitness=-838
## It 640: fitness=-838
## It 650: fitness=-838
## It 660: fitness=-838
## It 670: fitness=-838
## It 680: fitness=-838
## It 690: fitness=-838
## It 700: fitness=-838
## It 710: fitness=-838
## It 720: fitness=-838
## It 730: fitness=-838
## It 740: fitness=-838
## It 750: fitness=-838
## It 760: fitness=-838
## It 770: fitness=-838
## It 780: fitness=-838
## It 790: fitness=-838
## It 800: fitness=-838
## It 810: fitness=-838
## It 820: fitness=-838
## It 830: fitness=-838
## It 840: fitness=-838
## It 850: fitness=-838
## It 860: fitness=-838
## It 870: fitness=-838
## It 880: fitness=-838
## It 890: fitness=-838
## It 900: fitness=-838
## It 910: fitness=-838
## It 920: fitness=-838
## It 930: fitness=-838
## It 940: fitness=-838
## It 950: fitness=-838
## It 960: fitness=-838
## It 970: fitness=-838
## It 980: fitness=-838
## It 990: fitness=-838
## It 1000: fitness=-838
## Maximal number of iterations reached

Curvas de nivel

Optimizamos por gradiente multivariado

# Calculamos la derivada parcial
partial_dev <- function(x,i,fun,h=0.01){
  e <- x*0 # crea un vector de ceros de la misma longitud de x
  e[i] <- h
  y <- (fun(x+e)-fun(x-e))/(2*h)
  return(y)
}

# Obtenemos cada una de las derivadas parciales de f en x:
num_grad <- function(x,fun,h=0.01){
  # x: punto del espacio donde se debe evaluar el gradiente
  # fun: función para la que se desea calcular el gradiente en x
  # h: es el tamaño de ventana para el cálculo de la derivada numérica
  d <- length(x)
  y <- mapply(FUN=partial_dev,i=1:d,MoreArgs=list(x=x,h=h,fun=fun))
  return(y)
}

Sin precondicinamiento:

#pasamos schwefel al optimizador SIN precondicionamiento iniciando desde (380,390)
sol_schwefel <- optimizador_mult_numdev(f_schwefel2d_vec,x0=c(380,390),eta=1)

Curvas de nivel

contour(x1, x2, Z, las=1,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Función de Schwefel",
        sub = "Curvas de nivel de la función")
lines(sol_schwefel, type="b",cex=1.5,col="red")

Con preacondicionamiento:

#pasamos schwefel al optimizador CON precondicionamiento iniciando desde (380,390)
sol_schwefel_precon <- optimizador_mult_precond(f_schwefel2d_vec,x0=c(380,390),eta=1)

Curvas de nivel

contour(x1, x2, Z, las=1,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Función de Schwefel",
        sub = "Curvas de nivel de la función")
lines(sol_schwefel_precon, type="b",cex=1.5,col="red")
lines(sol_schwefel, type="b",cex=1.5,col="blue")
legend("topleft",col=c("red","blue"),legend = c("Con precondicionamiento","Sin precondicionamiento"),lty=1)

El problema del vendedor viajero

Para llevar acabo el desarrollo del problema del vendedor viajero primero tuvimos que crear la matriz de costos a partir de otras matrices. Obtuvimos la matriz “costos entre trayectos-peajes”, “kilometros entre ciudades” y “tiempo estimado entre ciudades en minutos”, como dos de estas matrices no estaban expresadas en terminos de costos, se multiplicaron por constantes como el sueldo por hora del vendedor, el precio de la gasolina por galon y el rendimiento de vehiculo, en este caso un Renault Sandero .

  • Sueldo basico del vendedor $10.417 pesos/hora
  • Precio del combustible $8.525 /Galon
  • Rendimiento del vehiculo 50 km/Galon

Dando como resultado la siguiente matriz de costos:

##              Palmira   Pasto   Tulua  Bogota Pereira Armenia Manizales
## Palmira        0.000 196.499  44.473 258.408 106.416  94.910   149.840
## Pasto        196.499   0.000 233.708 414.567 293.912 282.930   336.818
## Tulua         44.473 233.708   0.000 245.483  64.197  52.137   106.927
## Bogota       258.408 414.567 245.483   0.000 195.466 167.024   168.387
## Pereira      106.416 293.912  64.197 195.466   0.000  33.090    44.289
## Armenia       94.910 282.930  52.137 167.024  33.090   0.000    79.480
## Manizales    149.840 336.818 106.927 168.387  44.289  79.480     0.000
## Valledupar   518.590 682.718 489.496 385.731 415.293 509.128   371.515
## Monteria     469.152 637.597 425.718 445.872 365.858 408.135   360.014
## Soledad      574.551 737.117 542.333 465.389 503.189 550.466   460.435
## Cartagena    580.668 731.259 491.298 483.495 512.687 544.064   482.770
## Barranquilla 588.808 747.033 556.069 485.501 515.069 520.210   463.024
## Medellin     243.888 434.859 201.495 223.823 152.024 183.404   135.454
## Bucaramanga  365.393 519.523 315.651 203.123 266.435 295.972   234.164
## Cucuta       469.132 624.127 427.434 308.838 381.352 410.210   338.092
##              Valledupar Monteria Soledad Cartagena Barranquilla Medellin
## Palmira         518.590  469.152 574.551   580.668      588.808  243.888
## Pasto           682.718  637.597 737.117   731.259      747.033  434.859
## Tulua           489.496  425.718 542.333   491.298      556.069  201.495
## Bogota          385.731  445.872 465.389   483.495      485.501  223.823
## Pereira         415.293  365.858 503.189   512.687      515.069  152.024
## Armenia         509.128  408.135 550.466   544.064      520.210  183.404
## Manizales       371.515  360.014 460.435   482.770      463.024  135.454
## Valledupar        0.000  187.064 147.764   168.130      158.049  336.842
## Monteria        187.064    0.000 179.849   132.809      197.873  232.714
## Soledad         147.764  179.849   0.000    63.256        7.589  371.428
## Cartagena       168.130  132.809  63.256     0.000       72.780  350.760
## Barranquilla    158.049  197.873   7.589    72.780        0.000  387.881
## Medellin        336.842  232.714 371.428   350.760      387.881    0.000
## Bucaramanga     198.263  296.811 279.477   299.325      288.032  196.201
## Cucuta          220.555  312.853 300.563   327.173      315.707  303.415
##              Bucaramanga  Cucuta
## Palmira          365.393 469.132
## Pasto            519.523 624.127
## Tulua            315.651 427.434
## Bogota           203.123 308.838
## Pereira          266.435 381.352
## Armenia          295.972 410.210
## Manizales        234.164 338.092
## Valledupar       198.263 220.555
## Monteria         296.811 312.853
## Soledad          279.477 300.563
## Cartagena        299.325 327.173
## Barranquilla     288.032 315.707
## Medellin         196.201 303.415
## Bucaramanga        0.000 106.675
## Cucuta           106.675   0.000

Luego de obtener esta matriz procedemos a resolver el problema utilizando los dos metodos solicitados.

Algoritmos Geneticos (GA)

A continuacion presentamos el codigo utilizado para resolver el problema del vendedor viajero utilizando algoritmos geneticos, ademas incluye directamente el guardado de la grafica de cada resultado y genera un gif utilizando todas las graficas.

## Importamos las librerias que vamos a usar
library(GA)
library(tidyverse)
library(animation)
library(ggrepel)

# Creamos la matrix de coordenadas para las graficas
ciudades <- c("Palmira",    "Pasto",    "Tuluá",    "Bogota",   "Pereira",  "Armenia",  "Manizales",    "Valledupar",   "Montería", "Soledad",
              "Cartagena",  "Barranquilla", "Medellín", "Bucaramanga",  "Cúcuta")
coordenadas <- data.frame( long = c(-76.30361, -77.28111, -76.19536, -74.08175, -75.69611, -75.68111, -75.51738, -73.25322, -75.88143, -74.76459,
                                    -75.51444, -74.78132, -75.56359, -73.1198, -72.50782),
                           lat= c(3.53944, 1.21361, 4.08466, 4.60971, 4.81333, 4.53389, 5.06889, 10.46314, 8.74798, 10.91843,
                                  10.39972, 10.96854, 6.25184, 7.12539, 7.89391),
                           stringsAsFactors = F)

# Leemos la matrix de costos
D <- read.csv("costos_ciudades2.csv", sep = ";", row.name = 1)

#Funcion para calcular la longitud del tour
tourLength <- function(tour, distMatrix) {
  tour <- c(tour, tour[1])
  route <- embed(tour, 2)[,2:1]
  sum(distMatrix[route])
}

#Fitness function to be maximized
tspFitness <- function(tour, ...) 1/tourLength(tour, ...)

# Iteracciones maximas
Maxter <- c(25, 50, 75, 100, 200, 400, 600, 800, 1000,  1200)

# Ciudades
x <- coordenadas[,1]
y <- coordenadas[,2]

# Utilizamos la funcion saveGIF para crear el gif durante la iteraccion de la funcion de
# GA, es decir, durante cada entranamiento guardamos cada resultado en forma de grafica
# y por ultimos entregamos el GIF.
saveGIF(
  
  # For para cambiar el numero maximo de iteraciones en la funcion GA
  for (i in 1:length(Maxter)) {
    set.seed(124)

    GAiter <- ga(type = "permutation", fitness = tspFitness, distMatrix = D,
                 lower = 1, upper = 15, popSize = 50, maxiter = Maxter[i],
                 run = 500, pmutation = 0.2)

    tour <- GAiter@solution[1, ]
    tour <- c(tour, tour[1])
    n <- length(tour)

    ## Grafica
    map <- ggplot() +
      borders("world", "colombia") +
      geom_point(data = coordenadas,
                 aes(x = long,
                     y = lat), colour = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15))

    map_vectors <- map +
      ggtitle(paste( paste("Max. numero de iteraciones:", Maxter[i]), paste(" - Iteracciones realizadas: ", GAiter@iter))) +
      theme(plot.title = element_text(hjust = 0.5)) +
      geom_segment(aes(x = x[tour[-n]], y = y[tour[-n]], xend = x[tour[-1]], yend = y[tour[-1]]),
                   arrow = arrow(length = unit(0.25, "cm")), col = "steelblue", lwd= 0.6) +
      geom_text_repel(aes(x, y, label = ciudades),
                      size = 3.5)
    plot(map_vectors)

    # Guardamos las imagenes para el GIF
    ggsave(map_vectors, filename = gsub(" ", "", paste("Ga_", i,".png")),  width = 5.78 , height = 4.66)

  }
  , movie.name = "TSP_GA.gif")

A continuacion podremos obervar la primer grafica obtenida con el algoritmo genetico.

En esta grafica observamos que la ruta que siguio el algoritmo fue la siguiente:

##    PGA       PCitys
## 1   11    Cartagena
## 2   10      Soledad
## 3    9     Montería
## 4    5      Pereira
## 5    6      Armenia
## 6    3        Tuluá
## 7    1      Palmira
## 8    2        Pasto
## 9    7    Manizales
## 10   4       Bogota
## 11  13     Medellín
## 12  14  Bucaramanga
## 13   8   Valledupar
## 14  15       Cúcuta
## 15  12 Barranquilla

Dando como resultado un costo de:

## [1] 2667694

La siguiente es la ultima grafica donde podremos detallar que a pesar de tener un maximo de 1200 iteracion el algoritmo encontro la mejor ruta en tan solo 1044 iteraciones.

En esta grafica observamos que la ruta que siguio el algoritmo fue la siguiente:

##    BGA        Citys
## 1   11    Cartagena
## 2    9     Montería
## 3   13     Medellín
## 4    7    Manizales
## 5    5      Pereira
## 6    6      Armenia
## 7    3        Tuluá
## 8    1      Palmira
## 9    2        Pasto
## 10   4       Bogota
## 11  14  Bucaramanga
## 12  15       Cúcuta
## 13   8   Valledupar
## 14  10      Soledad
## 15  12 Barranquilla

Dando como resultado un costo de:

## [1] 2044516

Por ultimo, presentamos el gif en el cual observaremos las variaciones que el algoritmo presenta cuando varian el numero maximo de iteracciones posibles.

Colonia de Hormigas

Para este punto se llevo a acabo el desarrollo de la solucion en Python dado que para este lenguaje se presentaban una mayor cantidad de librerias disponibles, haciendo mas simple la implementacion del algoritmo y solucion del problema.

#Repo: https://github.com/johnberroa/Ant-Colony-Optimization
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import time
import warnings

from pandas.core.frame import DataFrame
warnings.filterwarnings("ignore")

class AntColonyOptimizer:
    def __init__(self, ants, evaporation_rate, intensification, alpha=1.0, beta=0.0, beta_evaporation_rate=0,
                 choose_best=.1):
        """
        Ant colony optimizer.  Traverses a graph and finds either the max or min distance between nodes.
        :param ants: number of ants to traverse the graph
        :param evaporation_rate: rate at which pheromone evaporates
        :param intensification: constant added to the best path
        :param alpha: weighting of pheromone
        :param beta: weighting of heuristic (1/distance)
        :param beta_evaporation_rate: rate at which beta decays (optional)
        :param choose_best: probability to choose the best route
        """
        # Parameters
        self.ants = ants
        self.evaporation_rate = evaporation_rate
        self.pheromone_intensification = intensification
        self.heuristic_alpha = alpha
        self.heuristic_beta = beta
        self.beta_evaporation_rate = beta_evaporation_rate
        self.choose_best = choose_best

        # Internal representations
        self.pheromone_matrix = None
        self.heuristic_matrix = None
        self.probability_matrix = None

        self.map = None
        self.set_of_available_nodes = None

        # Internal stats
        self.best_series = []
        self.best = None
        self.fitted = False
        self.best_path = None
        self.fit_time = None

        # Plotting values
        self.stopped_early = False

    def __str__(self):
        string = "Ant Colony Optimizer"
        string += "\n--------------------"
        string += "\nDesigned to optimize either the minimum or maximum distance between nodes in a square matrix that behaves like a distance matrix."
        string += "\n--------------------"
        string += "\nNumber of ants:\t\t\t\t{}".format(self.ants)
        string += "\nEvaporation rate:\t\t\t{}".format(self.evaporation_rate)
        string += "\nIntensification factor:\t\t{}".format(self.pheromone_intensification)
        string += "\nAlpha Heuristic:\t\t\t{}".format(self.heuristic_alpha)
        string += "\nBeta Heuristic:\t\t\t\t{}".format(self.heuristic_beta)
        string += "\nBeta Evaporation Rate:\t\t{}".format(self.beta_evaporation_rate)
        string += "\nChoose Best Percentage:\t\t{}".format(self.choose_best)
        string += "\n--------------------"
        string += "\nUSAGE:"
        string += "\nNumber of ants influences how many paths are explored each iteration."
        string += "\nThe alpha and beta heuristics affect how much influence the pheromones or the distance heuristic weigh an ants' decisions."
        string += "\nBeta evaporation reduces the influence of the heuristic over time."
        string += "\nChoose best is a percentage of how often an ant will choose the best route over probabilistically choosing a route based on pheromones."
        string += "\n--------------------"
        if self.fitted:
            string += "\n\nThis optimizer has been fitted."
        else:
            string += "\n\nThis optimizer has NOT been fitted."
        return string

    def _initialize(self):
        """
        Initializes the model by creating the various matrices and generating the list of available nodes
        """
        assert self.map.shape[0] == self.map.shape[1], "Map is not a distance matrix!"
        num_nodes = self.map.shape[0]
        self.pheromone_matrix = np.ones((num_nodes, num_nodes))
        # Remove the diagonal since there is no pheromone from node i to itself
        self.pheromone_matrix[np.eye(num_nodes) == 1] = 0
        self.heuristic_matrix = 1 / self.map
        self.probability_matrix = (self.pheromone_matrix ** self.heuristic_alpha) * (
                self.heuristic_matrix ** self.heuristic_beta)  # element by element multiplcation
        self.set_of_available_nodes = list(range(num_nodes))

    def _reinstate_nodes(self):
        """
        Resets available nodes to all nodes for the next iteration
        """
        self.set_of_available_nodes = list(range(self.map.shape[0]))

    def _update_probabilities(self):
        """
        After evaporation and intensification, the probability matrix needs to be updated.  This function
        does that.
        """
        self.probability_matrix = (self.pheromone_matrix ** self.heuristic_alpha) * (
                self.heuristic_matrix ** self.heuristic_beta)

    def _choose_next_node(self, from_node):
        """
        Chooses the next node based on probabilities.  If p < p_choose_best, then the best path is chosen, otherwise
        it is selected from a probability distribution weighted by the pheromone.
        :param from_node: the node the ant is coming from
        :return: index of the node the ant is going to
        """
        numerator = self.probability_matrix[from_node, self.set_of_available_nodes]
        if np.random.random() < self.choose_best:
            next_node = np.argmax(numerator)
        else:
            denominator = np.sum(numerator)
            probabilities = numerator / denominator
            next_node = np.random.choice(range(len(probabilities)), p=probabilities)
        return next_node

    def _remove_node(self, node):
        self.set_of_available_nodes.remove(node)

    def _evaluate(self, paths, mode):
        """
        Evaluates the solutions of the ants by adding up the distances between nodes.
        :param paths: solutions from the ants
        :param mode: max or min
        :return: x and y coordinates of the best path as a tuple, the best path, and the best score
        """
        scores = np.zeros(len(paths))
        coordinates_i = []
        coordinates_j = []
        for index, path in enumerate(paths):
            score = 0
            coords_i = []
            coords_j = []
            for i in range(len(path) - 1):
                coords_i.append(path[i])
                coords_j.append(path[i + 1])
                score += self.map[path[i], path[i + 1]]
            scores[index] = score
            coordinates_i.append(coords_i)
            coordinates_j.append(coords_j)
        if mode == 'min':
            best = np.argmin(scores)
        elif mode == 'max':
            best = np.argmax(scores)
        return (coordinates_i[best], coordinates_j[best]), paths[best], scores[best]

    def _evaporation(self):
        """
        Evaporate some pheromone as the inverse of the evaporation rate.  Also evaporates beta if desired.
        """
        self.pheromone_matrix *= (1 - self.evaporation_rate)
        self.heuristic_beta *= (1 - self.beta_evaporation_rate)

    def _intensify(self, best_coords):
        """
        Increases the pheromone by some scalar for the best route.
        :param best_coords: x and y (i and j) coordinates of the best route
        """
        i = best_coords[0]
        j = best_coords[1]
        self.pheromone_matrix[i, j] += self.pheromone_intensification

    def fit(self, map_matrix, iterations=100, mode='min', early_stopping_count=20, verbose=True):
        """
        Fits the ACO to a specific map.  This was designed with the Traveling Salesman problem in mind.
        :param map_matrix: Distance matrix or some other matrix with similar properties
        :param iterations: number of iterations
        :param mode: whether to get the minimum path or maximum path
        :param early_stopping_count: how many iterations of the same score to make the algorithm stop early
        :return: the best score
        """
        if verbose: print("Beginning ACO Optimization with {} iterations...".format(iterations))
        self.map = map_matrix
        start = time.time()
        self._initialize()
        num_equal = 0

        for i in range(iterations):
            start_iter = time.time()
            paths = []
            path = []

            for ant in range(self.ants):
                current_node = self.set_of_available_nodes[np.random.randint(0, len(self.set_of_available_nodes))]
                start_node = current_node
                while True:
                    path.append(current_node)
                    self._remove_node(current_node)
                    if len(self.set_of_available_nodes) != 0:
                        current_node_index = self._choose_next_node(current_node)
                        current_node = self.set_of_available_nodes[current_node_index]
                    else:
                        break

                path.append(start_node)  # go back to start
                self._reinstate_nodes()
                paths.append(path)
                path = []

            best_path_coords, best_path, best_score = self._evaluate(paths, mode)

            if i == 0:
                best_score_so_far = best_score
            else:
                if mode == 'min':
                    if best_score < best_score_so_far:
                        best_score_so_far = best_score
                        self.best_path = best_path
                elif mode == 'max':
                    if best_score > best_score_so_far:
                        best_score_so_far = best_score
                        self.best_path = best_path

            if best_score == best_score_so_far:
                num_equal += 1
            else:
                num_equal = 0

            self.best_series.append(best_score)
            self._evaporation()
            self._intensify(best_path_coords)
            self._update_probabilities()

            if verbose: print("Best score at iteration {}: {}; overall: {} ({}s)"
                              "".format(i, round(best_score, 2), round(best_score_so_far, 2),
                                        round(time.time() - start_iter)))

            if best_score == best_score_so_far and num_equal == early_stopping_count:
                self.stopped_early = True
                print("Stopping early due to {} iterations of the same score.".format(early_stopping_count))
                break

        self.fit_time = round(time.time() - start)
        self.fitted = True

        if mode == 'min':
            self.best = self.best_series[np.argmin(self.best_series)]
            if verbose: print(
                "ACO fitted.  Runtime: {} minutes.  Best score: {}".format(self.fit_time // 60, self.best))
            return self.best
        elif mode == 'max':
            self.best = self.best_series[np.argmax(self.best_series)]
            if verbose: print(
                "ACO fitted.  Runtime: {} minutes.  Best score: {}".format(self.fit_time // 60, self.best))
            return self.best
        else:
            raise ValueError("Invalid mode!  Choose 'min' or 'max'.")

    def plot(self):
        """
        Plots the score over time after the model has been fitted.
        :return: None if the model isn't fitted yet
        """
        if not self.fitted:
            print("Ant Colony Optimizer not fitted!  There exists nothing to plot.")
            return None
        else:
            fig, ax = plt.subplots(figsize=(20, 15))
            ax.plot(self.best_series, label="Best Run")
            ax.set_xlabel("Iteration")
            ax.set_ylabel("Performance")
            ax.text(.8, .6,
                    'Ants: {}\nEvap Rate: {}\nIntensify: {}\nAlpha: {}\nBeta: {}\nBeta Evap: {}\nChoose Best: {}\n\nFit Time: {}m{}'.format(
                        self.ants, self.evaporation_rate, self.pheromone_intensification, self.heuristic_alpha,
                        self.heuristic_beta, self.beta_evaporation_rate, self.choose_best, self.fit_time // 60,
                        ["\nStopped Early!" if self.stopped_early else ""][0]),
                    bbox={'facecolor': 'gray', 'alpha': 0.8, 'pad': 10}, transform=ax.transAxes)
            ax.legend()
            plt.title("Ant Colony Optimization Results (best: {})".format(np.round(self.best, 2)))
            plt.show()
#-------------------------------------------------------------------------------------------------------------

# Read the cost matrix
data = pd.read_csv('/home/joan/Escritorio/Ant_Colony_TSP/costos_ciudades.csv', sep=";", index_col=0)

# Change the dataframe to an array
data_matrix = data.to_numpy()

## CHANGE ANTS

# Create a dataframe to save iterations results
#best_pa = DataFrame(columns = ['Ants1','Ants15', 'Ants30', 'Ants45', 'Ants60', 'Ants75', 'Ants90','Ants105', 'Ants120', 'Ants135'])

# # Ants array                               
# Ants = [1, 15, 30, 45, 60, 75, 90, 105, 120, 135]

# for n in range(10):

#     # Create the Antiptimizer using differents parameteres
#     optimizer = AntColonyOptimizer(ants=Ants[n], evaporation_rate=.1, intensification=2, alpha=1, beta=1,
#                                    beta_evaporation_rate=0, choose_best=.1)

#     # Traing the model and save tha best result                             
#     optimizer.fit(data_matrix, 100, early_stopping_count = 20)
    
#     # Save the values in a dataframe
#     best_pa[best_pa.columns[n]] = optimizer.best_path

## CHANGE ITERATIONS

# Create a dataframe to save iterations results
#best_pa = DataFrame(columns = ['iter5','iter10', 'iter15', 'iter20', 'iter25', 'iter30', 'iter35', 'iter40', 'iter45', 'iter50'])
best_pa = DataFrame(columns = ['iter4','iter6', 'iter8', 'iter10', 'iter12', 'iter14', 'iter16', 'iter18', 'iter20', 'iter22'])

# Ants array                               
#Iter = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
Iter = [4, 6, 8, 10, 12, 14, 16, 18, 20, 22]

for n in range(10):

    # Create the Antiptimizer using differents parameteres
    optimizer = AntColonyOptimizer(ants=100, evaporation_rate=.1, intensification=2, alpha=1, beta=1,
                                   beta_evaporation_rate=0, choose_best=.1)

    # Traing the model and save tha best result                             
    optimizer.fit(data_matrix, Iter[n])
    
    # Save the values in a dataframe
    best_pa[best_pa.columns[n]] = optimizer.best_path

## SAVE RESULTS

# Show the best tour in terminal
print("\nTours table results: \n")
print(best_pa)

# saving the dataframe in a csv
print("\n Saving results in tours_results.csv")
best_pa.to_csv('/home/joan/Escritorio/Ant_Colony_TSP/tours_results.csv')

Para graficar y obtener el gif utilizamos el mismo codigo implementado con en algoritmos geneticos, solo que esta vez iteramos sobre el archivo de salida que genera el codigo en python.

## Usango colonia de hormigas

## Importamos las librerias que vamos a usar
#library(GA)
#library(tidyverse)
#library(animation)
#library(ggrepel)

# Creamos la matrix de coordenadas
#ciudades <- c("Palmira",   "Pasto",    "Tuluá",    "Bogota",   "Pereira",  "Armenia",  "Manizales",    "Valledupar",   "Montería", "Soledad",
#              "Cartagena", "Barranquilla", "Medellín", "Bucaramanga",  "Cúcuta")
#coordenadas <- data.frame( long = c(-76.30361, -77.28111, -76.19536, -74.08175, -75.69611, -75.68111, -75.51738, -73.25322, -75.88143, -74.76459,
#                                    -75.51444, -74.78132, -75.56359, -73.1198, -72.50782),
#                           lat= c(3.53944, 1.21361, 4.08466, 4.60971, 4.81333, 4.53389, 5.06889, 10.46314, 8.74798, 10.91843,
#                                  10.39972, 10.96854, 6.25184, 7.12539, 7.89391),
#                           stringsAsFactors = F)

# Leemos el archivo resultante luego de ejecutar el metodo de colonia de
# hormigas en python
Resultado_Ant <- read.csv("tours_results.csv", sep = ",", row.name = 1)
iter <- c(4,6,8,10,12,14,16,18,20,22)

# Ciudades
#x <- coordenadas[,1]
#y <- coordenadas[,2]

saveGIF(
  # For para cambiar los datos que vamos a graficar
  for (i in c(1,2,3,4,5,6,7,8,9,10)) {

    tour <- as.matrix(Resultado_Ant[i])
    tour <- tour + 1
    n <- length(tour)
    
    ## Grafica
    map <- ggplot() +
      borders("world", "colombia") +
      geom_point(data = coordenadas,
                 aes(x = long,
                     y = lat), colour = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15))
    
    map_vectors <- map +
      ggtitle(paste( paste("Numero de iteraciones:", iter[i]), "Hormigas: 100")) +
      theme(plot.title = element_text(hjust = 0.5)) +
      geom_segment(aes(x = x[tour[-n]], y = y[tour[-n]], xend = x[tour[-1]], yend = y[tour[-1]]),
                   arrow = arrow(length = unit(0.25, "cm")), col = "steelblue", lwd= 0.6) +
      geom_text_repel(aes(x, y, label = ciudades),
                      size = 3.5)
    plot(map_vectors)
    
    # Guardamos las imagenes para el GIF
    ggsave(map_vectors, filename = gsub(" ", "", paste("Ants_", i,".png")),  width = 5.78 , height = 4.66)
    
  }
  , movie.name = "TSP_ANTS.gif")

A continuacion podremos obervar la primer grafica obtenida con el algoritmo de colonia de hormigas

En esta grafica observamos que la ruta que siguio el algoritmo fue la siguiente:

##    PAN       PCitys
## 1    1      Palmira
## 2    3        Tuluá
## 3    6      Armenia
## 4    5      Pereira
## 5    7    Manizales
## 6   13     Medellín
## 7    4       Bogota
## 8   14  Bucaramanga
## 9   15       Cúcuta
## 10  10      Soledad
## 11  12 Barranquilla
## 12  11    Cartagena
## 13   8   Valledupar
## 14   9     Montería
## 15   2        Pasto

Dando como resultado un costo de:

## [1] 2413283

La siguiente es la ultima grafica donde podremos observar un parecido con la ultima grafica del algoritmo genetico.

En esta grafica observamos que la ruta que siguio el algoritmo fue la siguiente:

##    BAN        Citys
## 1   12 Barranquilla
## 2   10      Soledad
## 3    8   Valledupar
## 4   15       Cúcuta
## 5   14  Bucaramanga
## 6    4       Bogota
## 7    2        Pasto
## 8    1      Palmira
## 9    3        Tuluá
## 10   6      Armenia
## 11   5      Pereira
## 12   7    Manizales
## 13  13     Medellín
## 14   9     Montería
## 15  11    Cartagena

Dando como resultado un costo de:

## [1] 2044516

Por ultimo, presentamos el gif en el cual observaremos las variaciones que el algoritmo presenta cuando varian el numero maximo de iteracciones posibles.

Para concluir se puede observar que tanto en Algoritmos Geneticos como en Colonia de Hormigas el mejor costo obtenido para el problema del viajero vendedor fue de $2’044.516 pesos ademas que el recorrido de las ciudades es el mismo , lo que demuestra que ambos algoritmos convergen en la misma solucion y que son muy efectivos para la resolucion de este tipo de problemas.

Bibliografia