Parte 1: Optimización Numérica

Comenzamos importando las librerías necesarias para esta primera mitad.

library(lattice)
library(viridisLite)
library(GA)
library(pso)
library(DEoptim)
library(glue)
library(plotly)
library(ggplot2)
library(viridis)
library(reshape2)
library(animation)

Definición de funciones de gráficos y gifs para posterior uso simplificado.

# Función para simplificar producción gráfico 3D
plot3D_f <- function(x, y, z, input_title){
  # x, y : Coordenadas
  # z: Matrix con la función evaluada
  # input_title: Strin con el título
  plot_ly(
    x = x,
    y = y,
    z = z,
    type = "surface",
    # Escala de color moderna
    colorscale = "Viridis" 
  ) %>% layout(
    # Rotulación
    title = input_title,
    scene = list(
      xaxis = list(title = "X"),
      yaxis = list(title = "Y"),
      zaxis = list(title = "Z")
    )
  ) 
}

# Función para los gráficos de curvas de nivel
contourLevel <- function(x, y, z, nlevels, title, xlim = NULL, ylim = NULL){
  # x, y : Coordenadas
  # z: Matrix con la función evaluada
  # nlevels: Número de curvas de nivel a calcular
  # xlim, ylim: Dominio y Rango de la gráfica
  
  # Verificamos si se definieron límites
  if (!is.null(xlim) & !is.null(ylim)){
    contour(x, y, z, nlevels = nlevels, las = 1, xlim = xlim, ylim = ylim)
  }
  # Sino, se gráfica normalmente
  else{
    contour(x, y, z, nlevels = nlevels, las = 1)
  }
  # Rotulación
  title(main = title)
  mtext("X", side = 1, line = 3)
  mtext("Y", side = 2, line = 3, las = 1)
}

# Función para crear el GIF
makeGIF <- function(array_points, x_seq, y_seq, z_eval, num_levels,
                    title, file_name, interval, sol){
  # Definición de parámetros
    # array_points: Arreglo estructurado con los puntos de las poblaciones en cada iteración
    # x_seq: Secuencia de valores X
    # y_seq: Secuencia de valores Y
    # z_eval: Matriz con las evaluciones de la función para [x,y]
    # num_levels: Número de contornos de nivel
    # title: Título del gráfico
    # file_name: Nombre del archivo, con extensión .gif incluida
    # interval: Define los segundos entre frames
    # sol: Vector con las coordenadas de la solución real
  saveGIF({
    # Secuencia
    for (i in 1:length(array_points)) {
      # Frame
      contour(x_seq, y_seq, z_eval,
              nlevels = num_levels, col = "grey40",
              xlab = "X", ylab = "Y",
              main = paste(title, i))
      
      # Puntos
      points(array_points[[i]]$x, array_points[[i]]$y,
             col = "red", pch = 19)
      points(sol[1], sol[2],
             col = "green", pch = 19, lwd = 3)
    }
  # Renderización según el intervalo específicado y el nombre de archivo
  }, movie.name = file_name, interval = interval, ani.width = 400, ani.height = 400) 
}

Funciones de Prueba

Se seleccionaron las siguientes funciones de prueba para evaluar los diferentes algortimos de optimización:

Rastrigin.
Ronsenbrock.

La función de Rastrigin se define como: \[ f(\mathbf{x}) = An + \sum_{i=1}^{n} \left[ x_i^2 - A \cos(2\pi x_i) \right]\, \forall\,\, x_i \in [-5.12,\ 5.12] \] donde $A=10$ y $x_i \in [-5.12, 5.12]$. Dicha función cuenta con su mínimo global en $X = [0,0,...,0]$ con un valor de $f(\mathbf{x}) = 0$.

Luego, la de Rosenbrock se define como: \[ f(\mathbf{x}) = \sum_{i=1}^{n-1} \left[ 100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2 \right]\, \forall\,\, x_i \in [-\infty,\ \infty] \] Y esta función cuenta con su mínimo global en $X = [1,1,...,1]$ con un valor de $f(\mathbf{x}) = 0$.

Procedemos ahora a declarar las funciones vectorizadas en R.

# Función de Rastrigin con entrada como vector
rastrigin <- function(x) {
  A <- 10
  n <- length(x)
  return(A * n + sum(x^2 - A * cos(2 * pi * x)))
}
# Función de Rosenbrock con entrada como vector
rosenbrock <- function(x) {
  n <- length(x)
  # Aquí la indexación del vector de entrada
  # permite cumplir con la ecuación
  return(sum(100 * (x[-1] - x[-n]^2)^2 + (1 - x[-n])^2))
}
# Versión 1D -> 2D de Rosenbrock necesaria para algunos casos
rosenbrock_2d <- function(x) {
  (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
}

La función de Rastrigin tiene un dominio limitado, mientras la de Rosenbrock tiene como dominio los reales, sin embargo, ambas las vamos a trabajar en el espacio de búsqueda de Rastrigin, entendiendo que ambos mínimos globales están en dicho espacio, y ampliarlo es incrementar el costo computacional innecesariamente sin motivo alguno para el ejercicio de examinar los métodos de optimización.

# Definimos el espacio de búsqueda
x <- seq(-5.12, 5.12, length.out = 200)
y <- seq(-5.12, 5.12, length.out = 200)

# Cálculos 1D
f_ras <- sapply(x, function(val) rastrigin(c(val)))
f_ros <- sapply(x, function(val) rosenbrock_2d(c(val, 1)))

# Matrices para guardar valores de las funciones
z_ras <- matrix(0, length(x), length(y))
z_ros <- matrix(0, length(x), length(y))

# Calculamos los valores
for (i in 1:length(x)) {
  for (j in 1:length(y)) {
    # Ponemos la entrada como vector
    z_ras[i, j] <- rastrigin(c(x[i], y[j]))
    z_ros[i, j] <- rosenbrock(c(x[i], y[j]))
  }
}

# Creamos una grilla de coordenadas (x,y)
grid <- expand.grid(x = x, y = y)

# Evaluamos las funciones en toda la grilla de coordenadas
grid$z_ras <- apply(grid, 1, function(row) rastrigin(c(row["x"], row["y"])))
grid$z_ros <- apply(grid, 1, function(row) rosenbrock(c(row["x"], row["y"])))

Veamos ahora gráficamente como es la función de Rastrigin con entrada y salida 2D → 3D.

plot3D_f(x, y, z_ras, "Función de Rastrigin")

Fig.1 Superficie de Rastrigin

Igualmente un mapa de calor con curvas de nivel para entender desde otra perspectiva la función.

# Mapa de calor con contornos de Rastrigin
ggplot(grid, aes(x, y, z = z_ras)) +
  geom_raster(aes(fill = z_ras)) +
  geom_contour(color = "white", alpha = 0.5) +
  scale_fill_viridis_c() +
  labs(title = "Curvas de nivel - Rastrigin", x = "X", y = "Y", fill = "Z") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

Fig.2 Heatmap de Rastrigin

Ahora un gráfico de curvas de nivel simple para la función de Rastrigin.

contourLevel(x, y, z_ras, nlevels = 10, "Curvas de nivel de la función de Rastrigin")

Fig.3 Levelplot de Rastrigin

Hacemos lo mismo con la función de Rosenbrock para tener una idea general.

plot3D_f(x, y, z_ros, "Función de Rosenbrock")

Fig.4 Superficie de Rosenbrock

# Mapa de calor con contornos de Rosenbrock
ggplot(grid, aes(x, y, z = z_ros)) +
  geom_raster(aes(fill = z_ros)) +
  geom_contour(color = "white", alpha = 0.5) +
  scale_fill_viridis_c() +
  labs(title = "Curvas de nivel - Rosenbrock", x = "X", y = "Y", fill = "Z") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5))

Fig.5 Heatmap de Rosenbrock

contourLevel(x, y, z_ros, nlevels = 36, "Curvas de nivel de la función de Rosenbrock")

Fig.6 Levelplot de Rosenbrock

Con las funciones presentadas y teniendo claro su topografía, además del espacio de búsqueda; procedemos con las optimizaciones.

Optimización

NOTA: Para las funciones 3D → 4D no se hacen representaciones, debido a la abstracción fuera de nuestra percepción geométrica tridimensional.

Descenso por gradiente

El descenso de gradiente es un algoritmo de optimización determinista que busca minimizar (o maximizar) una función ajustando iterativamente sus variables en la dirección del gradiente negativo. Es ampliamente utilizado en matemáticas, ingeniería y aprendizaje automático como un método clásico. Tiene hiperparámetros como:

Punto inicial: Coordenadas desde donde comienza la búsqueda.
Tasa de aprendizaje (learning rate): Determina el tamaño del paso en cada iteración.
Número de iteraciones: Máximo de pasos antes de detenerse.
Criterio de convergencia: Umbral para detenerse si el cambio es muy pequeño.
Gradiente: Derivadas parciales de la función respecto a cada variable.

En esta parte se realizará la optimización de las anteriores funciones en dos y tres dimensiones, usando el método de descenso por gradiente, el cual consiste en recorrer la función dando pasos pequeños en la dirección en donde la función decrece más rápido.

En este caso debido a la complejidad de las funciones de prueba se realizaron algunas modificaciones al método simple, por ejemplo; si la distancia entre el punto y la actualización crece demasiado o el método no converge con cierto número de iteraciones, se reinicia el algoritmo con un nuevo punto aleatorio o el mejor punto encontrado hasta el momento, respectivamente. Esto hace que el método sea menos propenso a quedarse en mínimos locales y tenga más chances de llegar hasta el mínimo global.

Creamos entonces, la función para el decenso por gradiente

gradient_descent_mod <- function(f, grad_f, x0, eta = 0.001, tol = 1e-7, max_iter = 100000,
                                 max_restarts = 5, patience = 5000) {
  restart_count <- 0
  converged <- FALSE
  best_x <- x0
  best_value <- f(x0)
  last_improvement <- 0
  all_path <- list() # Lista para guardar los caminos de cada reinicio
  last_index <- 1

  while (!converged && restart_count < max_restarts) {
    x <- best_x # Reiniciar desde el mejor punto encontrado
    path <- matrix(NA, nrow = max_iter, ncol = length(x0))
    path[1, ] <- x # Guardar el punto inicial del reinicio

    cat("Reinicio número:", restart_count + 1, "Iniciando desde:", paste(x, collapse = ", "),
        "Valor:", best_value, "\n")

    for (i in 1:max_iter) {
      grad <- grad_f(x)
      last_index <- i + 1
      x_new <- x - eta * grad
      distance <- sqrt(sum((x_new - x)^2))

      if (is.infinite(distance)) {
        cat("La distancia es Inf. Reiniciando con un nuevo punto de inicio aleatorio.\n")
        x0_aleatorio <- runif(length(x0), min = -5, max = 5)
        best_x <- x0_aleatorio
        best_value <- f(x0_aleatorio)
        last_improvement <- i
        restart_count <- restart_count + 1
        break
      }

      path[i + 1, ] <- x_new # Guardar el nuevo punto
      value_new <- f(x_new)

      if (value_new < best_value) {
        best_value <- value_new
        best_x <- x_new
        last_improvement <- i
      }

      if (distance < tol) {
        path <- path[1:(i + 1), ]
        cat("Convergió en", i, "iteraciones (reinicio", restart_count, ")\n")
        converged <- TRUE
        break
      }

      # Criterio de paciencia para reiniciar
      if (i - last_improvement > patience && i > patience) {
        cat("No hay mejora significativa en", patience,
            "iteraciones. Reiniciando desde el mejor punto encontrado.\n")
        restart_count <- restart_count + 1
        break # Salir del bucle for y reiniciar
      }

      x <- x_new
    }

    all_path[[restart_count + 1]] <- path[1:last_index, ] # Guardar el camino del reinicio actual

    if (i == max_iter && !converged) {
      cat("No convergió en", max_iter, "iteraciones (reinicio", restart_count,
          "). Reiniciando desde el mejor punto encontrado.\n")
      restart_count <- restart_count + 1
    }
  }

  if (!converged) {
    cat("No se logró la convergencia después de", max_restarts, "reinicios.\n")
  }

  # Combinar todos los caminos en una sola matriz
  combined_path <- do.call(rbind, all_path)

  return(list(min = best_x, value = best_value, path = combined_path))
}

Rastrigin

Antes de proceder es importante definir la función gradiente de Rastrigin.

grad_rastrigin <- function(x) {
  A <- 10
  n <- length(x)
  grad <- numeric(n)
  for (i in 1:n) {
    grad[i] <- 2 * x[i] + 2 * pi * A * sin(2 * pi * x[i])
  }
  return(grad)
}

Comenzamos con la solución para 2 dimensiones base, o sea f(x,y):

# Punto inicial aleatorio
set.seed(2)
vec <- runif(2, min = -5.12, max = 5.12)
# Implementación del gradiente
gd_ras <- gradient_descent_mod(rastrigin, grad_rastrigin,
                                  x0 = vec, eta = 0.1, tol = 1e-3)

## Reinicio número: 1 Iniciando desde: -3.22680565834045, 2.07231012821198 Valor: 24.26904 
## No hay mejora significativa en 5000 iteraciones. Reiniciando desde el mejor punto encontrado.
## Reinicio número: 2 Iniciando desde: -1.02371988293758, 0.0575919115306958 Valor: 1.809776 
## No hay mejora significativa en 5000 iteraciones. Reiniciando desde el mejor punto encontrado.
## Reinicio número: 3 Iniciando desde: -1.01329607602703, -0.00972459521248048 Valor: 1.0804 
## No hay mejora significativa en 5000 iteraciones. Reiniciando desde el mejor punto encontrado.
## Reinicio número: 4 Iniciando desde: -1.01329607602703, -0.00972459521248048 Valor: 1.0804 
## No hay mejora significativa en 5000 iteraciones. Reiniciando desde el mejor punto encontrado.
## Reinicio número: 5 Iniciando desde: -1.01329607602703, -0.00972459521248048 Valor: 1.0804 
## No hay mejora significativa en 5000 iteraciones. Reiniciando desde el mejor punto encontrado.
## No se logró la convergencia después de 5 reinicios.

print(glue("\nLa mejor solución hallada está en (x, y) = ({gd_ras$min[1]}, {gd_ras$min[2]}).
Con un valor óptimo de {gd_ras$value}."))

## La mejor solución hallada está en (x, y) = (-1.01329607602703, -0.00972459521248048).
## Con un valor óptimo de 1.08040041990205.

Veamos como se ve la solución hallada en un gráfico de curvas de nivel, donde la solución óptima real está representada como un punto verde y la del método de optimización con un punto rojo:

contourLevel(x, y, z_ras, nlevels = 10, title = "Solución Descenso por Gradiente vs Exacta\nRastrigin",
             xlim = c(-3.5, 3.5), ylim = c(-3.5, 3.5))
points(gd_ras$min[1],gd_ras$min[2], pch = 19, col = "red", lwd = 1)
points(0,0, col = "green", pch = 19, lwd = 1)

Fig.7 GD Rastrigin

Se puede ver que a pesar de todos los reinicios e iteraciones del método, este se estancó en un mínimo local.

Veamos ahora para 3 dimensiones base, o sea f(x,y,z):

set.seed(2)
vec <- runif(3, min = -5.12, max = 5.12)
gd_ras <- gradient_descent_mod(rastrigin, grad_rastrigin,
                               x0 = c(3 , 3, 3), eta = 0.1,tol = 1e-3)

## Reinicio número: 1 Iniciando desde: 3, 3, 3 Valor: 27 
## No hay mejora significativa en 5000 iteraciones. Reiniciando desde el mejor punto encontrado.
## Reinicio número: 2 Iniciando desde: 4.58437385466581e-05, 4.58437385466581e-05, 4.58437385466581e-05 Valor: 1.250851e-06 
## Convergió en 6113 iteraciones (reinicio 1 )

print(glue("\nLa mejor solución hallada está en:
(x, y, z) = ({gd_ras$min[1]}, {gd_ras$min[2]}, {gd_ras$min[3]})
Con un valor óptimo de {gd_ras$value}"))

## La mejor solución hallada está en:
## (x, y, z) = (4.58437385466581e-05, 4.58437385466581e-05, 4.58437385466581e-05)
## Con un valor óptimo de 1.25085121283064e-06

Ese valor óptimo hallado en 3D, si analizamos las coordenadas (x, y, z) vemos que se aproxima al mínimo global, por lo que en esta ocasión el método lo hizo bien.

#### Rosenbrock

Nuevamente, antes de proceder con su optimización, formulamos la función gradiente de Rosenbrock.

# Gradiente de la función de Rosenbrock en 2D
grad_rosenbrock2D <- function(x) {
  dx <- -400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1])
  dy <- 200 * (x[2] - x[1]^2)
  return(c(dx, dy))
}

Aplicamos ahora si la optimización en 2D.

set.seed(2)
vec <- runif(2, min = -5.12, max = 5.12)
gd_ros <- gradient_descent_mod(rosenbrock, grad_rosenbrock2D, x0 = vec, eta = 0.001)

## Reinicio número: 1 Iniciando desde: -3.22680565834045, 2.07231012821198 Valor: 6973.367 
## La distancia es Inf. Reiniciando con un nuevo punto de inicio aleatorio.
## Reinicio número: 2 Iniciando desde: 0.733263348229229, -3.31948079634458 Valor: 1487.836 
## Convergió en 19788 iteraciones (reinicio 1 )

print(glue("\nLa mejor solución hallada está en (x, y) = ({gd_ros$min[1]}, {gd_ros$min[2]})
Con un valor óptimo de {gd_ros$value}"))

## La mejor solución hallada está en (x, y) = (0.99988826932461, 0.999776104001655)
## Con un valor óptimo de 1.25037364638396e-08

Veamos gráficamente el resultado para este caso, siguiendo la convención anterior para el óptimo real y el óptimo del método.

contourLevel(x, y, z_ros, nlevels = 4000, title = "Solución Descenso por Gradiente vs Exacta\nRosenbrock",
             xlim = c(0.5, 1.5), ylim = c(0.5, 1.5))
points(gd_ros$min[1],gd_ros$min[2], pch = 25, col = "red", lwd = 2)
points(1,1, col = "green", pch = 19, lwd = 1)

Fig.8 GD Rosenbrock

En este caso si convergió bien el método del gradiente con un punto inicial aletorio, posiblemente la topografía de la función de Rosenbrock no atrapa tan fácilmente con un mínimo local, o el punto inicial favoreció el proceso.

Ejecutemos ahora el caso 3D.

# Gradiente de la función de Rosenbrock en 3D
grad_rosenbrock3D <- function(x) {
  dx1 <- -400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1])
  dx2 <- 200 * (x[2] - x[1]^2) - 400 * x[2] * (x[3] - x[2]^2) - 2 * (1 - x[2])
  dx3 <- 200 * (x[3] - x[2]^2)
  return(c(dx1, dx2, dx3))
}

set.seed(2)
vec <- runif(3, min = -5.12, max = 5.12)
gd_ros <- gradient_descent_mod(rosenbrock, grad_rosenbrock3D, x0 = vec, eta = 0.001)

## Reinicio número: 1 Iniciando desde: -3.22680565834045, 2.07231012821198, 0.750861668586731 Valor: 8230.232 
## La distancia es Inf. Reiniciando con un nuevo punto de inicio aleatorio.
## Reinicio número: 2 Iniciando desde: -3.31948079634458, 4.43839338840917, 4.43474958650768 Valor: 27661.62 
## La distancia es Inf. Reiniciando con un nuevo punto de inicio aleatorio.
## Reinicio número: 3 Iniciando desde: -3.70841023279354, 3.33448815625161, -0.319814844988286 Valor: 23964.93 
## La distancia es Inf. Reiniciando con un nuevo punto de inicio aleatorio.
## Reinicio número: 4 Iniciando desde: 0.499837417155504, 0.526740669738501, -2.6110524055548 Valor: 842.4896 
## Convergió en 17095 iteraciones (reinicio 3 )

print(glue("\nLa mejor solución hallada está en:
(x, y, z) = ({gd_ros$min[1]}, {gd_ros$min[2]}, {gd_ros$min[3]})
Con un valor óptimo de {gd_ros$value}"))

## La mejor solución hallada está en:
## (x, y, z) = (0.999954281647075, 0.999908391106169, 0.999816354226281)
## Con un valor óptimo de 1.05044371078583e-08

Según los resultamos numéricos, nuevamente convergió bien el método.

NOTA: Para los siguientes métodos heurísticos que se basan en poblaciones e iteraciones, estos parámetros serán iguales para todos; 50 miembros y 100 iteraciones.

Algoritmos evolutivos

Los algoritmos evolutivos (AE) son métodos de optimización inspirados en los mecanismos de la evolución biológica. Funcionan con una población de posibles soluciones que evolucionan a lo largo de varias generaciones mediante operaciones como la selección, cruzamiento, mutación y elitismo. La función ga de la librería GA la implementa y cuenta con parámetros como:

Población inicial (popSize): Número de individuos (soluciones candidatas) generadas al inicio.
Función de aptitud (fitness): Evalúa qué tan buena es una solución. Si el problema es de minimización, se suele usar el valor negativo de la función.
Cruce o recombinación (pcrossover): Combina dos soluciones para crear nuevas.
Mutación (pmutation): Modifica aleatoriamente una solución para mantener diversidad.
Elitismo (elitism): Número de mejores soluciones que se copian directamente a la siguiente generación.
Generaciones máximas (maxiter): Número total de generaciones.
Criterio de parada temprana (run): Número de generaciones sin mejora antes de detenerse.
Optimización local opcional (optim): Puede afinar soluciones localmente.

El objetivo es encontrar una solución óptima o cercana al óptimo para un problema, especialmente cuando los métodos tradicionales no son eficaces debido a la complejidad del espacio de búsqueda.

Rastrigin

Comenzamos con el caso 1D.

genetic_rastri <- ga(
  type = "real-valued",
  fitness = function(x){return(-rastrigin(x))},
  lower = c(-5.12),
  upper = c(5.12),
  # Población
  popSize = 50,
  # Iteraciones máximas permitidas
  maxiter = 100,
  optim = TRUE,
  # Semilla generadora de pseudo-aleatorios
  seed = 2,
  monitor = FALSE
)

# Mostrar la mejor solución y su fitness:
print(glue("Mejor solución encontrada:
x = {genetic_rastri@solution[1]}
Con un valor óptimo de {genetic_rastri@fitnessValue} en {genetic_rastri@iter} iteraciones"))

## Mejor solución encontrada:
## x = 0
## Con un valor óptimo de 0 en 100 iteraciones

# Graficar la curva
plot(x, f_ras, type = "l", col = "blue", lwd = 2,
     main = "Solución Algoritmo Evolutivo vs Exacta\nRastrigin",
     xlab = "x", ylab = "f(x)")
grid()
points(genetic_rastri@solution[1], genetic_rastri@fitnessValue, col = "red", pch = 25, lwd = 3)
points(0, 0, col = "green", pch = 19, lwd = 1)

Fig.9 GA Rastrigin 1D

Vemos que no supone mayor dificultad para el método hallar una solución exacta en baja dimensionalidad.

Procedemos con la de Rastrigin 2D, considerando diversos argumentos de la función.

genetic_rastri <- ga(
  type = "real-valued",
  fitness = function(x){return(-rastrigin(x))},
  lower = c(-5.12, -5.12),
  upper = c(5.12, 5.12),
  # Población
  popSize = 50,
  # Iteraciones máximas permitidas
  maxiter = 100,
  optim = TRUE,
  # Semilla generadora de pseudo-aleatorios
  seed = 2,
  monitor = FALSE
)

# Mostrar la mejor solución y su fitness:
print(glue("Mejor solución encontrada:
(x, y) = ({genetic_rastri@solution[1]}, {genetic_rastri@solution[2]})
Con un valor óptimo de {genetic_rastri@fitnessValue} en {genetic_rastri@iter} iteraciones"))

## Mejor solución encontrada:
## (x, y) = (0, 0)
## Con un valor óptimo de 0 en 100 iteraciones

contourLevel(x, y, z_ras, nlevels = 30, "Solución Algoritmo Evolutivo vs Exacta\nRastrigin",
             xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(genetic_rastri@solution[1],genetic_rastri@solution[2],
       pch = 25, col = "red", lwd = 2)
points(0, 0, col = "green", pch = 19, lwd = 1)

Fig.10 GA Rastrigin 2D

Nuevamente, la solución hallada es exacta para este caso de 2 dimensiones de entrada.

Veamos que tal con la versión 3D.

genetic_rastri <- ga(
  type = "real-valued",
  fitness = function(x) -rastrigin(x),
  lower = c(-5.12, -5.12, -5.12),
  upper = c(5.12, 5.12, 5.12),
  popSize = 50,
  maxiter = 100,
  optim = TRUE,
  seed = 2,
  monitor = FALSE
)

# Mostrar la mejor solución y su fitness:
print(glue("Mejor solución encontrada:
(x, y, z) = ({genetic_rastri@solution[1]}, {genetic_rastri@solution[2]}, {genetic_rastri@solution[3]})
Con un valor óptimo de {genetic_rastri@fitnessValue} en {genetic_rastri@iter} iteraciones"))

## Mejor solución encontrada:
## (x, y, z) = (0, 0, 0)
## Con un valor óptimo de 0 en 100 iteraciones

De nuevo, solución perfectamente exacta.

Rosenbrock

Es momento de evaluar ahora que tal le va al método con la de Rosenbrock, si logrará la misma exactitud que con Rastrigin.

genetic_rosen <- ga(
  type = "real-valued",
  fitness = function(x){return(-rosenbrock_2d(c(x, 1)))},
  lower = c(-5.12),
  upper = c(5.12),
  # Población
  popSize = 50,
  # Iteraciones máximas permitidas
  maxiter = 100,
  optim = TRUE,
  # Semilla generadora de pseudo-aleatorios
  seed = 2,
  monitor = FALSE
)

# Mostrar la mejor solución y su fitness:
print(glue("Mejor solución encontrada:
x = {genetic_rosen@solution[1]}
Con un valor óptimo de {genetic_rosen@fitnessValue} en {genetic_rosen@iter} iteraciones"))

## Mejor solución encontrada:
## x = 0.999999501247124
## Con un valor óptimo de -9.97504772334171e-11 en 100 iteraciones

# Graficar la curva
plot(x, f_ros, type = "l", col = "blue", lwd = 2,
     main = "Solución Algoritmo Evolutivo vs Exacta\nRosenbrock",
     xlab = "x", ylab = "f(x)", ylim = c(0, 10))
grid()
points(genetic_rosen@solution[1], genetic_rosen@fitnessValue, col = "red", pch = 25, lwd = 3)
points(1, 0, col = "green", pch = 19, lwd = 1)

Fig.11 GA Rosenbrock 1D

Se obtiene una solución que a términos prácticos es exacta.

Luego, para el caso 2D de Rosenbrock

genetic_rosen <- ga(
  type = "real-valued",
  fitness = function(x) -rosenbrock(x),
  lower = c(-5.12, -5.12),
  upper = c(5.12, 5.12),
  popSize = 50,
  maxiter = 100,
  optim = TRUE,
  seed = 2,
  monitor = FALSE
)

# Mostrar la mejor solución y su fitness:
print(glue("Mejor solución encontrada:
(x, y) = ({genetic_rosen@solution[1]}, {genetic_rosen@solution[2]})
Con un valor óptimo de {genetic_rosen@fitnessValue} en {genetic_rosen@iter} iteraciones"))

## Mejor solución encontrada:
## (x, y) = (0.999821444661359, 0.9996430029726)
## Con un valor óptimo de -3.18826775556571e-08 en 100 iteraciones

contourLevel(x, y, z_ros, nlevels = 4000, "Solución Algoritmo Evolutivo vs Exacta\nRosenbrock",
             xlim = c(0.5, 1.5), ylim = c(0.5, 1.5))
points(genetic_rosen@solution[1],genetic_rosen@solution[2],
       pch = 25, col = "red", lwd = 2)
points(1, 1, col = "green", pch = 19, lwd = 1)

Fig.12 GA Rosenbrock 2D

Continua siendo muy buena la solución sin problema alguno.

Ahora el de Rosenbrock 3D.

genetic_rosen <- ga(
  type = "real-valued",
  fitness = function(x) -rosenbrock(x),
  lower = c(-5.12, -5.12, -5.12),
  upper = c(5.12, 5.12, 5.12),
  popSize = 50,
  maxiter = 100,
  optim = TRUE,
  seed = 2,
  monitor = FALSE
)

# Mostrar la mejor solución y su fitness:
print(glue("Mejor solución encontrada:
(x, y, z) = ({genetic_rosen@solution[1]}, {genetic_rosen@solution[2]}, {genetic_rosen@solution[3]})
Con un valor óptimo de {genetic_rosen@fitnessValue} en {genetic_rosen@iter} iteraciones"))

## Mejor solución encontrada:
## (x, y, z) = (0.999902313108414, 0.999805319634613, 0.999610896342258)
## Con un valor óptimo de -4.74947456265742e-08 en 100 iteraciones

Según los resultados númericos, se aproxima bastante bien.

Enjambres de partículas

Se basa en las ideas de comportamientos de poblaciones como grupos de aves, que en búsqueda de alimento, no solo se mueven a cierta velocidad, sino que buscan en base a lo que observan del entorno, y lo que observan en sus compañeros. Para aplicar la optimización por enjambre de partículas se empleó la función psoptim de la librería pso (particle swarm optimization). Esta permite modificar diversos parámetros como por ejemplo:

Máximo número de iteraciones.
Tamaño del enjambre.
Magnitud del vector de velocidad.
Constante de exploración local (cognitivo).
Constante de exploración global (social).

Entre muchos otros parámetros adicionales, pero, para este caso particular de comparación de métodos de optimización, solo se modificará el máximo número de iteraciones y población, los demás parámetros se dejarán por defecto.

Rastrigin

set.seed(2)
# Usamos la función con 100 iteraciones
ras_optim = psoptim(par = rep(NA, 2), fn = rastrigin,
                    control = list(maxit = 100, s = 50),
                    lower = c(-5.12), upper = c(5.12))

print(glue("La mejor solución hallada es:
x = {ras_optim$par[1]}, con f(x) = {ras_optim$value}.
Usando {ras_optim$counts['iteration']} iteraciones."))

## La mejor solución hallada es:
## x = 1.73824685829077e-05, con f(x) = 6.00024101515828e-08.
## Usando 100 iteraciones.

# Graficar la curva
plot(x, f_ras, type = "l", col = "blue", lwd = 2,
     main = "Solución PSO vs Exacta\nRastrigin",
     xlab = "x", ylab = "f(x)")
grid()
points(ras_optim$par[1], ras_optim$value, col = "red", pch = 25, lwd = 3)
points(0, 0, col = "green", pch = 19, lwd = 1)

Fig.13 PSO Rastrigin 1D

La solucion númerica se aproxima lo suficientemente bien para ser gráficamente perfecta.

set.seed(2)
# Usamos la función con 100 iteraciones y 50 miembros
ras_optim = psoptim(par = rep(NA, 2), fn = rastrigin,
                    control = list(maxit = 100, s = 50),
                    lower = c(-5.12, -5.12), upper = c(5.12, 5.12))

print(glue("La mejor solución hallada es:
x = {ras_optim$par[1]}, y = {ras_optim$par[2]}, con f(x,y) = {ras_optim$value}.
Usando {ras_optim$counts['iteration']} iteraciones."))

## La mejor solución hallada es:
## x = 1.73824685829077e-05, y = 5.41619145005576e-07, con f(x,y) = 6.00024101515828e-08.
## Usando 100 iteraciones.

Veamos en una gráfica de curvas de nivel la solución hallada:

contourLevel(x, y, z_ras, nlevels = 30, "Solución PSO vs Exacta\nRastrigin",
             xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(ras_optim$par[1], ras_optim$par[2], pch = 25, col = "red", lwd = 2)
points(0,0, col = "green", pch = 19, lwd = 1)

Fig.14 PSO Rastrigin 2D

Continuamos con una solución óptima del método, suficientemente buena gráficamente, con un error despreciable.

set.seed(2)
# Usamos la función con 100 iteraciones
ras_optim = psoptim(par = rep(NA, 3), fn = rastrigin,
                    control = list(maxit = 100, s = 50),
                    lower = rep(-5.12, 3), upper = rep(5.12, 3))

print(glue("La mejor solución hallada es:
           x = {ras_optim$par[1]}, y = {ras_optim$par[2]}, con f(x,y) = {ras_optim$value}.
           Usando {ras_optim$counts['iteration']} iteraciones."))

## La mejor solución hallada es:
## x = 0.00795172284832642, y = -0.000250197075416108, con f(x,y) = 0.0126238210993428.
## Usando 100 iteraciones.

Se sigue aproximando decentemente, pero ya empieza a resultarle más complicado al método.

Rosenbrock

set.seed(2)
# Usamos la función con 100 iteraciones
ros_optim = psoptim(par = rep(NA, 2), fn = rosenbrock_2d,
                    control = list(maxit = 100, s = 50),
                    lower = c(-5.12), upper = c(5.12))

print(glue("La mejor solución hallada es:
           x = {ras_optim$par[1]}, con f(x) = {ras_optim$value}.
           Usando {ras_optim$counts['iteration']} iteraciones."))

## La mejor solución hallada es:
## x = 0.00795172284832642, con f(x) = 0.0126238210993428.
## Usando 100 iteraciones.

# Graficar la curva
plot(x, f_ros, type = "l", col = "blue", lwd = 2,
     main = "Solución PSO vs Exacta\nRosenbrock",
     xlab = "x", ylab = "f(x)", ylim = c(0, 10))
grid()
points(ros_optim$par[1], ros_optim$value, col = "red", pch = 25, lwd = 3)
points(1, 0, col = "green", pch = 19, lwd = 1)

Fig.15 PSO Rosenbrock 1D

No le cuesta mucho en 1D aproximarse bien, pero de primeras no logra el mismo nivel de exactitud numérico que los algoritmos evolutivos.

Ahora veamos el resultado de optimizar la de Rosenbrock 2 variables.

set.seed(2)
# Usamos la función con 100 iteraciones
ros_optim = psoptim(par = rep(NA, 2), fn = rosenbrock,
                    control = list(maxit = 100, s = 50),
                    lower = rep(-5.12, 2), upper = rep(5.12, 2))

print(glue("La mejor solución hallada es:
x = {ros_optim$par[1]}, y = {ros_optim$par[2]}, con f(x,y) = {ros_optim$value}.
Usando {ras_optim$counts['iteration']} iteraciones."))

## La mejor solución hallada es:
## x = 1.00096511939299, y = 1.0018622281904, con f(x,y) = 1.40675608269388e-06.
## Usando 100 iteraciones.

contourLevel(x, y, z_ros, nlevels = 4000, "Solución PSO vs Exacta\nRosenbrock",
             xlim = c(0.5, 1.5), ylim = c(0.5, 1.5))
points(ros_optim$par[1], ros_optim$par[2], pch = 25, col = "red", lwd = 2)
points(1,1, col = "green", pch = 19, lwd = 1)

Fig.16 PSO Rosenbrock 2D

Ahora veamos el resultado de optimizar la de Rosenbrock con variables x,y,z.

set.seed(2)
# Usamos la función con 100 iteraciones
ros_optim = psoptim(par = rep(NA, 3), fn = rosenbrock,
                    control = list(maxit = 100, s = 50),
                    lower = rep(-5.12, 3), upper = rep(5.12, 3))

print(glue("La mejor solución hallada es:
           x = {ros_optim$par[1]}, y = {ros_optim$par[2]}, z = {ros_optim$par[3]}, con f(x,y,z) = {ros_optim$value}.
           Usando {ras_optim$counts['iteration']} iteraciones."))

## La mejor solución hallada es:
## x = 0.898746387671411, y = 0.814071484869055, z = 0.66124977362888, con f(x,y,z) = 0.0490379825415431.
## Usando 100 iteraciones.

Vemos que ya en 3D el método de enjambre de partículas, le cuesta obtener resultados decentes con la función de Rosenbrock, se evidencia principalmente en la 3ra coordenada, con $z\approx 0.6612$.

Evolución diferencial

La evolución diferencial (DE) es un algoritmo evolutivo específico, que se diferencia por su metodología para generar variación entre individuos por medio de diferencias entre ellos, modulada por un factor de escala F, a parte de tener las otras operaciones de los algortimos evolutivos convencionales. Los parámetros principales de este algoritmo incluyen:

Tamaño de población.
Factor de escala.
Probabilidad de cruce.
Número de generaciones.

DE es adecuado para problemas con espacios de búsqueda no lineales y multidimensionales, y es menos propenso a quedar atrapado en óptimos locales en comparación con otros algoritmos de optimización.

Comenzamos creando un objeto de control del algoritmo, ya que la librería lo requiere de esa forma previamente.

# Creamos el objeto de control necesario para la función DEoptim
ctrl <- DEoptim.control(
  NP      = 50,    # tamaño de población
  itermax = 100,   # iteraciones máximas
  CR      = 0.9,   # tasa de cruce
  F       = 0.8,   # factor de escala
  trace   = FALSE  # no imprimir cada generación
)

Rastrigin

Iniciamos esta sub-sección con Rastrigin 1D.

set.seed(2)
de_rastri <- DEoptim(
  fn      = rastrigin,
  lower   = c(-5.12),
  upper   = c(5.12),
  control = ctrl
)

print(glue("Valor óptimo Rastrigin 1D → f* = {de_rastri$optim$bestval}.
           Obtenido en x = {de_rastri$optim$bestmem[1]}.
           Usando {de_rastri$optim$iter} iteraciones."))

## Valor óptimo Rastrigin 1D → f* = 1.04805053524615e-13.
## Obtenido en x = 2.2989377927243e-08.
## Usando 100 iteraciones.

# Graficar la curva
plot(x, f_ras, type = "l", col = "blue", lwd = 2,
     main = "Solución Evolución Diferencial vs Exacta\nRastrigin",
     xlab = "x", ylab = "f(x)")
grid()
points(de_rastri$optim$bestmem[1], de_rastri$optim$bestval, col = "red", pch = 25, lwd = 3)
points(0, 0, col = "green", pch = 19, lwd = 1)

Fig.17 DE Rastrigin 1D

Empieza bien, pero no precisamente con la mayor de las exactitudes de todos los métodos.

Proseguimos optimizando la función en 2D.

set.seed(2)
de_rastri <- DEoptim(
  fn      = rastrigin,
  lower   = c(-5.12, -5.12),
  upper   = c(5.12, 5.12),
  control = ctrl
)

print(glue("Valor óptimo Rastrigin 2D → f* = {de_rastri$optim$bestval}.
Obtenido en (x,y) = ({de_rastri$optim$bestmem[1]}, {de_rastri$optim$bestmem[2]}).
Usando {de_rastri$optim$iter} iteraciones."))

## Valor óptimo Rastrigin 2D → f* = 0.00154272911834141.
## Obtenido en (x,y) = (0.00172512919681201, 0.00219093500219472).
## Usando 100 iteraciones.

Evaluemos la calidad de la solución visualmente con una gráfica de contornos de nivel.

best_point <- de_rastri$optim$bestmem

contourLevel(x, y, z_ras, nlevels = 30, "Solución Evolución Diferencial vs Exacta\nRastrigin",
             xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(best_point[1], best_point[2], pch = 25, col = "red", lwd = 2)
points(0, 0, col = "green", pch = 19, lwd = 1)

Fig.18 DE Rastrigin 2D

Aquí ya se aprecia ligeramente mayor desviación en la solución hallada.

Veamos ahora como le va para el caso de tres variables de entrada.

# Optimización de la función Rastrigin en 3D
set.seed(2)
de_rastri<- DEoptim(
  fn = rastrigin,
  lower = rep(-5.12, 3),
  upper = rep(5.12, 3),
  control = ctrl
)

print(glue("Valor óptimo Rastrigin 3D → f* = {de_rastri$optim$bestval}.
Obtenido en (x,y,z) = ({de_rastri$optim$bestmem[1]}, {de_rastri$optim$bestmem[2]}, {de_rastri$optim$bestmem[3]}).
Usando {de_rastri$optim$iter} iteraciones."))

## Valor óptimo Rastrigin 3D → f* = 0.191570036875426.
## Obtenido en (x,y,z) = (0.0223378538065587, 0.0214201641130212, 0.00305181153790221).
## Usando 100 iteraciones.

Como era de esperar, cada vez le resulta más costoso mejorar la solución en Rastrigin al aumentar la dimensionalidad.

Rosenborck

Continuamos ahora con la optimización por evolución diferencial para la de Rosenbrock en 1D.

set.seed(2)
de_rosen <- DEoptim(
  fn      = function(x){return(rosenbrock_2d(c(x,1)))},
  lower   = c(-5.12),
  upper   = c(5.12),
  control = ctrl
)

print(glue("Valor óptimo Rosenbrock 1D → f* = {de_rosen$optim$bestval}.
Obtenido en x = {de_rosen$optim$bestmem[1]}.
Usando {de_rosen$optim$iter} iteraciones."))

## Valor óptimo Rosenbrock 1D → f* = 3.34126966787017e-27.
## Obtenido en x = 1.
## Usando 100 iteraciones.

# Graficar la curva
plot(x, f_ros, type = "l", col = "blue", lwd = 2,
     main = "Solución Evolución Diferencial vs Exacta\nRosenbrock",
     xlab = "x", ylab = "f(x)", ylim = c(0, 10))
grid()
points(de_rosen$optim$bestmem[1], de_rosen$optim$bestval, col = "red", pch = 25, lwd = 3)
points(1, 0, col = "green", pch = 19, lwd = 1)

Fig.19 DE Rosenbrock 1D

Vemos que la evolución diferencial arranca teniendo un resultado exacto en el argumento que minimiza a la función de Rosenbrock.

Ahora bien, veamos si el buen desempeño se mantiene con una entrada 2D.

set.seed(2)
de_rosen <- DEoptim(
  fn      = rosenbrock,
  lower   = c(-5.12, -5.12),
  upper   = c(5.12, 5.12),
  control = ctrl
)

print(glue("Valor óptimo Rosenbrock 2D → f* = {de_rosen$optim$bestval}.
Obtenido en (x,y) = ({de_rosen$optim$bestmem[1]}, {de_rosen$optim$bestmem[2]}).
Usando {de_rosen$optim$iter} iteraciones."))

## Valor óptimo Rosenbrock 2D → f* = 9.87537269674006e-21.
## Obtenido en (x,y) = (1.00000000006936, 1.00000000014584).
## Usando 100 iteraciones.

best_point <- de_rosen$optim$bestmem

contourLevel(x, y, z_ros, nlevels = 4000, "Solución Evolución Diferencial vs Exacta\nRosenbrock",
              xlim = c(0.5, 1.5), ylim = c(0.5, 1.5))
points(best_point[1], best_point[2], pch = 25, col = "red", lwd = 2)
points(1, 1, col = "green", pch = 19, lwd = 1)

Fig.20 DE Rosenbrock 2D

El método continua produciendo un muy buen resultado con esta función.

Veamos que tal con entrada 3D (x, y, z).

# Optimización de la función Rosenbrock en 3D
set.seed(2)
de_rosen <- DEoptim(
  fn = rosenbrock,
  lower = rep(-5.12, 3),
  upper = rep(5.12, 3),
  control = ctrl
)
print(glue("Valor óptimo Rosenbrock 3D → f* = {de_rosen$optim$bestval}.
Obtenido en (x,y,z) = ({de_rosen$optim$bestmem[1]}, {de_rosen$optim$bestmem[2]}, {de_rosen$optim$bestmem[3]}).
Usando {de_rosen$optim$iter} iteraciones."))

## Valor óptimo Rosenbrock 3D → f* = 2.00531391570777e-11.
## Obtenido en (x,y,z) = (0.99999829525011, 0.999996667066797, 0.999993567646815).
## Usando 100 iteraciones.

A términos generales, la evolución diferencial le va especialmente bien con la función de Rosenbrock

Generación de los GIF

Para poder generar los gifs, es necesario activar el monitoreo o trazabilidad de los miembros de las poblaciones, para cada iteración, por lo cual se volverán a correr todos los métodos de optmización pero con parámetros específicos.

Observación: Las celdas de producción de GIFs solo se ejecutaron una vez para generar los archivos, y ya luego se etiquetaron con HTML. Para la ejecución del .Rmd cambiar “eval = FALSE” a “eval = TRUE”, o en su defecto borrar ese parámetro del chunk.

Descenso por gradiente

Comenzamos produciendo el GIF para Rastrigin.

set.seed(2)
vec <- runif(2, min = -5.12, max = 5.12)
gd_ras_gif <- gradient_descent_mod(rastrigin, grad_rastrigin, x0 = vec, eta = 0.001)

positions <- list()
for (iter in 1:dim(gd_ras_gif$path)[1]) {
  df <- data.frame(x = gd_ras_gif$path[iter,1], y = gd_ras_gif$path[iter,2])
  # Guardo los puntos de esa iteración en la lista
  positions[[iter]] <- df
}

# Generamos el GIF
makeGIF(positions, x, y, z_ras, 10, "Descenso por Gradiente en Rastrigin - Iteración",
        "GD_rastrigin2D.gif", interval = 0.2, sol = c(0, 0))

Gif.1 Proceso Gradiente Rastrigin

Se aprecia inmediatamente que el método convergió muy rápido, pero a un mínimo local.

Miramos el GIF de Rosenbrock con gradiente.

set.seed(2)
vec <- runif(2, min = -5.12, max = 5.12)
gd_ros_gif <- gradient_descent_mod(rosenbrock, grad_rosenbrock2D, x0 = vec, eta = 0.001)

positions <- list()
for (iter in 1:dim(gd_ros_gif$path)[1]) {
  df <- data.frame(x = gd_ros_gif$path[iter,1], y = gd_ros_gif$path[iter,2])
  # Guardo los puntos de esa iteración en la lista
  positions[[iter]] <- df
}

# Antes de pasar al GIF, como este proceso tuvo demasiadas iteraciones
# y el método de gradiente se estanca mucho, tomamos puntos cada 200 muestras
positions <- positions[seq(1, length(positions), by = 200)]

# Generamos el GIF
num_iter <- dim(gd_ros_gif$path)[1]
makeGIF(positions, x, y, z_ros, 40,
        "Descenso por Gradiente en Rosenbrock - Iteración",
        "GD_rosenbrock2D.gif", interval = 0.2, sol = c(1, 1))

Gif.2 Proceso Gradiente Rosenbrock

Recordando que acorde al código, debido a la enorme cantidad de iteraciones del método de gradiente, se tomaron iteraciones en pasos grandes; se puede observar que aunque son muchas iteraciones, logra alcanzar el mínimo global. Esto sumado a las ejecuciones iniciales de experimentación, denota que la función de Rosenbrock, tiene una topografía que no perjudica tanto al método, si lo comparamos con la de Rastrigin por ejemplo.

**NOTA: Para los GIF de los métodos heurísticos, se emplean 20 miembros de población y 100 iteraciones.*

Algoritmos evolutivos

Para poder guardar los puntos de cada población en cada iteracioń con la librería GA es necesario crear previamente una función de monitoreo, y luego pasarla a los argumentos de la función ga.

# Creamos lista global para guardar los puntos por generación
poblacion_por_gen <- list()

# Monitor personalizado que guarda la población tras cada generación
mi_monitor <- function(obj) {
  poblacion_por_gen[[obj@iter + 1]] <<- obj@population  # Guardar población
}

# Ejecutar algoritmo genético
genetic_rastrigin_gif <- ga(
  type = "real-valued",
  fitness = function(x){return(-rastrigin(x))},   # Minimización
  lower = c(-5.12, -5.12),
  upper = c(5.12, 5.12),
  popSize = 20,
  maxiter = 100,
  run = 100,
  monitor = mi_monitor
)

# Guardamos los puntos pero omitimos el primer elemento que siempre es NULL en GA
poblacion_por_gen <- poblacion_por_gen[2:101]
positions <- list()
for (iter in 1:length(poblacion_por_gen)) {
  mat <- t(poblacion_por_gen[[iter]])
  df <- data.frame(x = mat[1, ], y = mat[2, ])
  positions[[iter]] <- df
}

# Generamos el gif
makeGIF(positions, x, y, z_ras, 10, "GA en Rastrigin - Iteración",
        "GA_rastrigin2D.gif", interval = 0.2, sol = c(0, 0))

Gif.3 Proceso Evolutivo Rastrigin

Se aprecia que hay puntos que llegan al mínimo global, pero al final de las iteraciones, varios miembros permanecieron en mínimos locales, aquí la ventaja de obtener la buena solución es por el almacenamiento de las coordenadas que dan el mejor valor detectado al final del algoritmo.

Veamos ahora para Rosenbrock con GA.

# Creamos lista global para guardar los puntos por generación
poblacion_por_gen <- list()

# Monitor personalizado que guarda la población tras cada generación
mi_monitor <- function(obj) {
  poblacion_por_gen[[obj@iter + 1]] <<- obj@population  # Guardar población
}

# Ejecutar algoritmo genético
genetic_rosenbrock_gif <- ga(
  type = "real-valued",
  fitness = function(x) -rosenbrock(x),   # Minimización
  lower = c(-5.12, -5.12),
  upper = c(5.12, 5.12),
  popSize = 20,
  maxiter = 100,
  run = 100,
  monitor = mi_monitor
)

# Guardamos los puntos pero omitimos el primer elemento que siempre es NULL en GA
poblacion_por_gen <- poblacion_por_gen[2:101]
positions <- list()
for (iter in 1:length(poblacion_por_gen)) {
  mat <- t(poblacion_por_gen[[iter]])
  df <- data.frame(x = mat[1, ], y = mat[2, ])
  positions[[iter]] <- df
}

# Generamos el gif
makeGIF(positions, x, y, z_ros, 40, "GA en Rosenbrock - Iteración",
        "GA_rosenbrock2D.gif", interval = 0.2, sol = c(1, 1))

Gif.4 Proceso Evolutivo Rosenbrock

Se aprecia que los miembros de la población a medida que pasan las iteraciones, se van acercando muy suavemente al óptimo, sin dejar de existir puntos que realizan exploraciones para no sesgar la población en un camino evolutivo.

Enjambres de partículas

Para poder activar el monitoreo que trae la función de librería, es necesario setear varios argumentos,los cuales son:

trace = 1.
REPORT = 1.
trace_stats = TRUE.

Si alguno de estos no se tiene activado, el monitoreo solo tomará datos de cada cierto número de iteraciones, o tomará solo los valores de las soluciones y no los puntos.

set.seed(2)
result_rastrigin <- psoptim(
  par = rep(NA, 2),
  fn = rastrigin,
  lower = rep(-5.12, 2),
  upper = rep(5.12, 2),
  control = list(
    maxit = 100,
    s = 20,
    trace = 1,
    REPORT = 1,
    trace.stats = TRUE
  )
)

# Guardar los puntos de las partículas
positions <- list()
for (iter in 1:length(result_rastrigin$stats$x)) {
  mat <- result_rastrigin$stats$x[[iter]]
  df <- data.frame(x = mat[1, ], y = mat[2, ])
  positions[[iter]] <- df
}

makeGIF(positions, x, y, z_ras, 10, "PSO en Rastrigin - Iteración",
        "PSO_rastrigin2D.gif", interval = 0.2, sol = c(0, 0))

Gif.5 Proceso Enjambre Rastrigin

El método converge bastante bien, sin dejar de tener individuos en búsqueda constante.

Produzcamos ahora el GIF para Rosenbrock.

set.seed(2)
result_rosenbrock <- psoptim(
  par = rep(NA, 2),
  fn = rosenbrock,
  lower = rep(-5.12, 2),
  upper = rep(5.12, 2),
  control = list(
    maxit = 100,
    s = 20,
    trace = 1,
    REPORT = 1,
    trace.stats = TRUE
  )
)
positions <- list()
for (iter in 1:length(result_rosenbrock$stats$x)) {
  mat <- result_rosenbrock$stats$x[[iter]]
  df <- data.frame(x = mat[1, ], y = mat[2, ])
  positions[[iter]] <- df
}

makeGIF(positions, x, y, z_ros, 40, "PSO en Rosenbrock - Iteración",
        "PSO_rosenbrock2D.gif", interval = 0.2, sol = c(1, 1))

Gif.6 Proceso Enjambre Rosenbrock

Particularmente en la función de Rosenbrock, aunque el método de enjambre va convergiendo, varios miembros se acumulan alrededor del óptimo, por cuestiones de la naturaleza cognitiva de los miemrbos del algoritmo.

Evolución diferencial

Esta función solo requiere setear los argumentos storepopfrom y storepopfreq para guardar los puntos de población en cada iteración, y se accede a ellos por medio del atributo $member del objeto.

Comenzamos primero, con la estructura de control.

ctrl <- DEoptim.control(
  NP      = 20,
  itermax = 100,
  CR      = 0.9,
  F       = 0.8, 
  trace   = FALSE,  
  storepopfrom = 1,
  storepopfreq = 1,
)

set.seed(2) 
de_rastri <- DEoptim(
  fn      = rastrigin,
  lower   = c(-5.12, -5.12),
  upper   = c(5.12, 5.12),
  control = ctrl
)

pop_list <- de_rastri$member$storepop

positions <- vector("list", length(pop_list))
for (i in seq_along(pop_list)) {
  pop_mat <- pop_list[[i]]
  positions[[i]] <- data.frame(x = pop_mat[,1],
                                        y = pop_mat[,2])
}

makeGIF(positions, x, y, z_ras, 10, "DE en Rastrigin - Iteración",
        "DE_rastrigin2D.gif", interval = 0.2, sol = c(0, 0))

Gif.7 Proceso Diferencial Rastrigin

Por ejemplo en este GIF se aprecia la importancia del tamaño de población, sin individuos suficientes en constante búsqueda, el algoritmo se estancaría fácilmente en un óptimo local. Aquí también se aprecia como miembros de la población van migrando a diversos mínimos de la función, haciendo uso del distintivo del algoritmo, que es la comparación por diferencias.

Finalizamos los GIFs con el de Rosenbrock.

set.seed(2) 
de_rosen <- DEoptim(
  fn      = rosenbrock,
  lower   = c(-5.12, -5.12),
  upper   = c(5.12, 5.12),
  control = ctrl
)

pop_list <- de_rosen$member$storepop

positions <- vector("list", length(pop_list))
for (i in seq_along(pop_list)) {
  pop_mat <- pop_list[[i]]
  positions[[i]] <- data.frame(x = pop_mat[,1],
                                        y = pop_mat[,2])
}

makeGIF(positions, x, y, z_ros, 40, "DE en Rosenbrock - Iteración",
        "DE_rosenbrock2D.gif", interval = 0.2, sol = c(1, 1))

Gif.8 Proceso Diferencial Rosenbrock

Se resalta la alta, rápida y organizada convergencia del método de evolución diferencial con la función de Rosenbrock, bajo los hiperparámetros estándar establecidos.

Heuristícos VS Descenso por gradiente

Para tener una herramienta previo al análisis final, se va establecer una métrica comparativa (entendiendo las simplificaciones que esto implica por la diversidad de métodos), en la que se va a medir el error en optimización promedio de cada método, para cada función, y este error se va a ponderar en función de la cantidad media de evaluaciones de la F.O e iteraciones del método. Más explicitamente:

\[ \text{Métrica de ineficiencia} = \overline{|{\text{error}}|} \times \overline{\text{iteraciones}} \times \overline{\text{evaluaciones}} \] Para esto vamos a correr en un ciclo definido for, los 3 métodos heurísticos con 100 semillas diferentes de aleatoriedad, recopilaremos los resultados, obtendremos el promedio de cada factor, y con ellos calcularemos la métrica. La prueba será con las funciones a 3 variables de entrada, para mayor exigencia, y como parámetros lo estableceremos con 50 miembros de población y 100 iteraciones límite.

Rastrigin

n_runs <- 100

# Estructura de control para evolución diferencial
ctrl <- DEoptim.control(
  NP      = 50,
  itermax = 100,
  trace   = FALSE,  
)

# Listas para los datos
ga_error <- vector("list", 100)
ga_iter <- vector("list", 100)
ga_eval <- vector("list", 100)

pso_error <- vector("list", 100)
pso_iter <- vector("list", 100)
pso_eval <- vector("list", 100)

de_error <- vector("list", 100)
de_iter <- vector("list", 100)
de_eval <- vector("list", 100)

# Simulación con 100 semillas diferentes
for (i in 1:n_runs) {
  # GA
  set.seed(i)
  ga_i <- ga(
    type = "real-valued",
    fitness = function(x){return(-rastrigin(x))},
    lower = rep(-5.12, 3),
    upper = rep(5.12, 3),
    popSize = 50,
    maxiter = 100,
    monitor = FALSE
  )
  ga_error[[i]] <- abs(ga_i@fitnessValue)
  ga_iter[[i]] <- ga_i@iter
  # No hay método directo, toca calcularlo así
  ga_eval[[i]] <- ga_i@popSize * ga_i@iter
  
  # PSO
  set.seed(i)
  pso_i <- psoptim(
    par = rep(NA, 3),
    fn = rastrigin,
    control = list(maxit = 100, s = 50),
    lower = rep(-5.12, 3), upper = rep(5.12, 3)
  )
  pso_error[[i]] <- abs(pso_i$value)
  pso_iter[[i]] <- pso_i$counts["iteration"]
  pso_eval[[i]] <- pso_i$counts["function"]
  
  # DE
  set.seed(i)
  de_i <- DEoptim(
    fn      = rastrigin,
    lower   = rep(-5.12, 3), upper = rep(5.12, 3),
    control = ctrl
  )
  de_error[[i]] <- abs(de_i$optim$bestval)
  de_iter[[i]] <- de_i$optim$iter
  de_eval[[i]] <- de_i$optim$nfeval
}

# Calculamos la métrica
metric_ga = prod(sapply(list(ga_error, ga_iter, ga_eval), function(x) mean(unlist(x))))
metric_pso = prod(sapply(list(pso_error, pso_iter, pso_eval), function(x) mean(unlist(x))))
metric_de = prod(sapply(list(de_error, de_iter, de_eval), function(x) mean(unlist(x))))

Veamos gráficamente los resultados de la simulación con Rastrigin.

# Creamos un data frame
datos <- data.frame(
  categoria = c("GA", "PSO", "DE"),
  valor = round(c(metric_ga, metric_pso, metric_de),3)
)

colors <- c("GA" = "#A9DFBF", "PSO" = "#AED6F1", "DE" = "#FFD7B5")

# Gráfico de barras
ggplot(datos, aes(x = categoria, y = valor, fill = categoria)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = valor), vjust = -0.5) +
  scale_fill_manual(values = colors) +
  labs(title = "Métrica de ineficiencia - Rastrigin",
       x = "Métodos",
       y = "Ineficiencia") +
  theme_minimal() +
  theme(legend.position = "none")

Fig.21 Comparación Rastrigin

Rosenbrock

# Listas para los datos
ga_error <- vector("list", 100)
ga_iter <- vector("list", 100)
ga_eval <- vector("list", 100)

pso_error <- vector("list", 100)
pso_iter <- vector("list", 100)
pso_eval <- vector("list", 100)

de_error <- vector("list", 100)
de_iter <- vector("list", 100)
de_eval <- vector("list", 100)

# Simulación con 100 semillas diferentes
for (i in 1:n_runs) {
  # GA
  set.seed(i)
  ga_i <- ga(
    type = "real-valued",
    fitness = function(x){return(-rosenbrock(x))},
    lower = rep(-5.12, 3),
    upper = rep(5.12, 3),
    popSize = 50,
    maxiter = 100,
    monitor = FALSE
  )
  ga_error[[i]] <- abs(ga_i@fitnessValue)
  ga_iter[[i]] <- ga_i@iter
  # No hay método directo, toca calcularlo así
  ga_eval[[i]] <- ga_i@popSize * ga_i@iter
  
  # PSO
  set.seed(i)
  pso_i <- psoptim(
    par = rep(NA, 3),
    fn = rosenbrock,
    control = list(maxit = 100, s = 50),
    lower = rep(-5.12, 3), upper = rep(5.12, 3)
  )
  pso_error[[i]] <- abs(pso_i$value)
  pso_iter[[i]] <- pso_i$counts["iteration"]
  pso_eval[[i]] <- pso_i$counts["function"]
  
  # DE
  set.seed(i)
  de_i <- DEoptim(
    fn      = rosenbrock,
    lower   = rep(-5.12, 3), upper = rep(5.12, 3),
    control = ctrl
  )
  de_error[[i]] <- abs(de_i$optim$bestval)
  de_iter[[i]] <- de_i$optim$iter
  de_eval[[i]] <- de_i$optim$nfeval
}

# Calculamos la métrica
metric_ga = prod(sapply(list(ga_error, ga_iter, ga_eval), function(x) mean(unlist(x))))
metric_pso = prod(sapply(list(pso_error, pso_iter, pso_eval), function(x) mean(unlist(x))))
metric_de = prod(sapply(list(de_error, de_iter, de_eval), function(x) mean(unlist(x))))

Veamos gráficamente los resultados de la iteración con Rosenbrock.

# Creamos un data frame
datos <- data.frame(
  categoria = c("GA", "PSO", "DE"),
  valor = round(log(c(metric_ga, metric_pso, metric_de),3))
)

colors <- c("GA" = "#A9DFBF", "PSO" = "#AED6F1", "DE" = "#FFD7B5")

# Gráfico de barras
ggplot(datos, aes(x = categoria, y = valor, fill = categoria)) +
  geom_bar(stat = "identity") +
  geom_text(aes(label = valor), vjust = -0.5) +
  scale_fill_manual(values = colors) +
  labs(title = "Métrica de ineficiencia - Rosenbrock (Logarítmica)",
       x = "Métodos",
       y = "Ineficiencia") +
  theme_minimal() +
  theme(legend.position = "none")

Fig.22 Comparación Rosenbrock

Fue necesario poner esta gráfica con el logaritmo base euler de los valores, ya que la métrica de ineficiencia de algoritmos evolutivos GA era mucho mayor a las demás.

Discusión

Tras las ejecuciones realizadas tanto en materia de experimentaciones iniciales, ejecuciones oficiales en el documento, como los análisis a partir de los GIF; todo esto en complemento con la teoría conocida del curso y las discuiones de clase, se puede concluir lo siguiente:

El método de optimización por descenso de gradiente es muy dependiente de sus hiperparámetros y de la topografía del cuerpo matemático a optimizar, siendo los argumentos más relevantes:
- Punto inicial: Es clave, pues un punto inicial muy bueno o conveniente puede significar la convergencia al óptimo global y un menor costo computacional por iteraciones.
- Learning Rate: Según la naturaleza del problema, puede implicar alejarse del óptimo global buscado y dificultar su convergencia debido a “movimientos bruscos” en el espacio de búsqueda. También puede implicar un descenso demasiado lento, lo cual conlleva más iteraciones y un alto riesgo de estancarse en un óptimo local.
El punto inicial en el descenso por gradiente puede tener dos facetas. En ejercicios didácticos o de revisión, elegir un punto conveniente puede percibirse como hacer trampa. En cambio, en ejercicios complejos con fines de investigación, puede basarse en un conocimiento profundo sobre la naturaleza del problema, permitiendo seleccionar un punto inicial con criterios robustos.
Los métodos heurísticos empleados se basan en ideas bioinspiradas poderosas. Esto se entiende al considerar que la naturaleza, con el paso del tiempo, ha llevado a los seres vivos a optimizar sus procesos con fines de supervivencia y reproducción. En ese sentido, los algoritmos de este tipo tienen una base sólida respaldada por la propia naturaleza, a pesar de las simplificaciones presentes en algunas de sus implementaciones.
Una diferencia clave que explica el éxito de los métodos heurísticos por encima del descenso por gradiente es la implementación de poblaciones, ya que estas ofrecen diferentes frentes de búsqueda de óptimos, minimizando el riesgo de quedar atrapados en óptimos locales o de requerir iteraciones extensivas. Los miembros de la población se respaldan y guían entre sí.
Los algoritmos bioinspirados demuestran ser una estrategia computacional versátil y potente para abordar problemas de optimización complejos, especialmente cuando la información del gradiente no está disponible o el espacio de búsqueda exhibe una naturaleza extensa y no lineal que dificulta obtener una expresión analítica.
La fortaleza puntual de los algoritmos evolutivos está en la exploración de sus individuos potenciada por la combinación de sus atributos, mediante cruce y la introducción de diversidad a través de mutación.
Equiparar o hacer comparaciones extensas entre los propios métodos heurísticos representa una tarea complicada, ya que poseen hiperparámetros muy distintos entre sí, los cuales no se pueden igualar de forma simple. Ni siquiera la cantidad de iteraciones máximas puede compararse directamente, dado que una iteración en un algoritmo no tiene el mismo peso o significado que en otro. Además, de la dependencia del estado aleatorio de inicialización de la población en menor medida.
A pesar de la dificultad para equipararlos de forma detallada; con las simulaciones realizadas se pudo evidenciar, bajo algunos casos estándar, ciertas virtudes de unos métodos sobre otros:
- Por ejemplo, la evolución diferencial logra resultados eficientes con la función de Rosenbrock, pues fue la de menor métrica de ineficiencia, mientras que a los algoritmos evolutivos les costó dramáticamente dicha función.
- Por otro lado y para contraste, los algoritmos evolutivos fueron los mejores para Rastrigin, mientras que los enjambres de partículas resultaron ser los más ineficientes para esta función. Esto deja a la evolución diferencial para estas simulaciones estándar, como una opción altamente equilibrada.

Parte 2: Optimización Combinatoria

En esta sección, se aborda la optimización del recorrido de un vendedor que debe visitar 13 ciudades principales de Colombia: Bogotá, Cali, Medellín, Barranquilla, Cartagena, Cúcuta, Bucaramanga, Pereira, Santa Marta, Ibagué, Pasto, Manizales y Neiva. El objetivo central es determinar la ruta que minimice el costo total del viaje, empleando dos metaheurísticas: Optimización por Colonias de Hormigas (ACO) y Algoritmos Genéticos (GA). El costo total considera factores como el salario del conductor durante el viaje, el consumo de combustible, las distancias, tiempos de viaje y peajes.

Descripción del problema del vendedor

El problema se modela como una instancia del Problema del Viajante de Comercio (TSP), un clásico desafío de optimización combinatoria. Se busca encontrar la ruta cíclica más económica que visite cada una de las 13 ciudades exactamente una vez y regrese a la ciudad de origen. El costo asociado al desplazamiento entre cualquier par de ciudades se desglosa en tres componentes principales:

Costo del salario : Se tomó como referencia el salario mínimo mensual legal vigente (SMMLV) en Colombia al momento del estudio, fijado en $1,423,500 COP (fuente: Ministerio de Trabajo). Haciendo cálculos simples llegamos a que la hora ordinaria en Colombia se da por la siguiente fórmula:

\[ \text{Valor hora ordinaria} = \frac{\text{salario mensual}}{\left(\frac{\text{horas semanales}}{6} \times 30\right)} \]

Las horas semanales en Colombia a fecha de la realización de este trabajo es de 46, reemplazando en la ecuación:

\[ \text{Valor hora ordinaria} = \frac{$1.423.500}{\left(\frac{46}{6} \times 30\right)} \]

Lo cual nos da un valor de hora ordinaria de $6.189, el cual usaremos para el calculo del costo del salario.

Adicionalmente necesitamos el tiempo de conducción entre ciudades para poder calcular el costo del salario. Para esto, se consultaron en Google Maps Distance Matrix API utilizando el modo driving.

El costo del salario se calcula con la siguiente ecuación (1):

$\text{Costo_salario} = \text{tiempo} \times \text{costo_hora} \tag{1}$
Costo del combustible . Se utilizó un precio promedio del galón de gasolina, obtenido a partir de los precios reportados por la Comisión de Regulación de Energía y Gas (CREG) para las 13 ciudades principales en la fecha del análisis, resultando en $15,827 COP/galón.

Se seleccionó como vehículo de referencia un Kia Picanto 1.0L. Según datos de la publicación especializada “El Carro Colombiano”, este vehículo presenta un consumo mixto aproximado de 60 km/galón.

Además necesitamos la distancia entre las ciudades, la cual fue obtenida de la misma API de Google Maps. La distancia entre las ciudades se obtiene en kilómetros y se usa para calcular el costo del combustible.

El costo del combustible se calcula con la siguiente ecuación (2):

\[ \text{Costo_gasolina} = \left( \frac{\text{Distancia}}{\text{Rendimiento (km/galon)}} \right) \times \text{Precio_por_galon} \tag{2} \]

Observación: Se uso el galon como unidad de medida ya que en Colombia es el estandar para medir la gasolina.
Costo de peajes. La información sobre los costos de peajes para las rutas que conectan las 13 ciudades se obtuvo de la plataforma Tollguru (Se consideraron los peajes de las rutas más comunes).

El costo de los peajes se calcula con la siguiente ecuación (3):

\[ \text{Costo_peaje} = \text{peaje} \tag{3} \]

La suma de todos estos costos nos da el costo total de desplazamiento entre dos ciudades, que se calcula con la siguiente ecuación (4): \[ \text{Costo_total} = \text{Costo_salario} + \text{Costo_gasolina} + \text{Costo_peaje} \tag{4} \] Este problema se modela como un Problema del Viajante de Comercio (TSP, por sus siglas en inglés), en el cual buscamos la ruta más corta o económica para visitar todas las ciudades.

Algoritmos utilizados

En este caso, se usan dos algoritmos de optimización para encontrar la mejor ruta:

1. Algoritmo de Colonias de Hormigas (ACO): Este algoritmo se inspira en el comportamiento colectivo de las hormigas al buscar caminos. Emplea feromonas artificiales para marcar las rutas: las rutas más cortas (de menor costo) acumulan más feromonas, incrementando la probabilidad de que sean elegidas por hormigas futuras, guiando así la búsqueda hacia soluciones prometedoras.

2. Algoritmo Genético (GA): Basado en los principios de la evolución biológica y la selección natural, este algoritmo opera sobre una población de posibles rutas (cromosomas). Mediante operadores como la selección (favoreciendo las rutas de menor costo), el cruce (combinando partes de rutas existentes para crear nuevas) y la mutación (introduciendo pequeñas alteraciones aleatorias), la población evoluciona a lo largo de generaciones, convergiendo hacia soluciones de alta calidad.

Código en R para la optimización

Cargamos las librerias necesarias:

library(stats)
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sf)

## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE

library(rnaturalearth)
library(rnaturalearthdata)

## 
## Attaching package: 'rnaturalearthdata'

## The following object is masked from 'package:rnaturalearth':
## 
##     countries110

library(gganimate)
library(gifski)
library(GA)
library(httr)

## 
## Attaching package: 'httr'

## The following object is masked from 'package:plotly':
## 
##     config

library(jsonlite)
library(dplyr)
library(tibble)
library(readr)
library(dotenv)

Consultamos las distancias viales entre las ciudades y tiempos de viaje ofrecidos por la API de Google. Adicionalmente cargamos los costos de peajes, los nombres de las ciudades y el número de ciudades:

# Consultamos las distancias viales y tiempos de viaje
dotenv::load_dot_env()
API_KEY <- Sys.getenv("APIKEY_DE_GOOGLE_MAPS")

addresses <- c(
  "Bogota, Colombia", "Cali, Colombia", "Medellin, Colombia", "Barranquilla, Colombia",
  "Cartagena, Colombia", "Cucuta, Colombia", "Bucaramanga, Colombia", "Pereira, Colombia",
  "Santa Marta, Colombia", "Ibague, Colombia", "Pasto, Colombia", "Manizales, Colombia",
  "Neiva, Colombia"
)

city_names <- sub(",.*$", "", addresses)


get_distance_time_matrix <- function(addr_vec, labels, api_key) {
  n <- length(addr_vec)
  dist_mat <- matrix(NA_real_, n, n)
  time_mat <- matrix(NA_real_, n, n)
  
  for (i in seq_along(addr_vec)) {
    res <- httr::GET(
      "https://maps.googleapis.com/maps/api/distancematrix/json",
      query = list(
        origins      = addr_vec[i],
        destinations = paste(addr_vec, collapse = "|"),
        key          = api_key,
        mode         = "driving",
        language     = "es",
        units        = "metric"
      )
    )
    httr::stop_for_status(res)
    
    # 👇  ¡OJO!  simplifyVector = FALSE mantiene la jerarquía original
    json <- fromJSON(content(res, "text", encoding = "UTF-8"),
                     simplifyVector = FALSE)
    
    if (json$status != "OK") {
      stop("API error: ", json$error_message)
    }
    
    elements <- json$rows[[1]]$elements   # ahora sí es lista de 13 sub-listas
    
    # Recorremos destino por destino
    for (j in seq_along(elements)) {
      el <- elements[[j]]
      if (el$status == "OK") {
        dist_mat[i, j] <- el$distance$value / 1000
        time_mat[i, j] <- el$duration$value / 3600
      }
    }
    
    Sys.sleep(1)
  }
  
  dist_df <- as_tibble(dist_mat, .name_repair = ~labels) |> 
             add_column(Ciudad = labels, .before = 1)
  time_df <- as_tibble(time_mat, .name_repair = ~labels) |> 
             add_column(Ciudad = labels, .before = 1)
  dist_df <- as.data.frame(dist_mat)
  time_df <- as.data.frame(time_mat)

  rownames(dist_df) <- labels
  rownames(time_df) <- labels
  colnames(dist_df) <- labels
  colnames(time_df) <- labels
  list(distances = dist_df, times = time_df)
}

matrices <- get_distance_time_matrix(addresses, city_names, API_KEY)
df_distancias  <- matrices$distances
df_tiempos_viajes  <- matrices$times

# Guardar a CSV
write_csv(df_distancias,  "df_distancias.csv")
write_csv(df_tiempos_viajes,  "df_tiempos_viajes.csv")


df_peajes <- read.csv("df_peajes.csv", row.names = 1)
ciudades_nombres <- c("Bogota", "Cali", "Medellin", "Barranquilla", "Cartagena", 
                      "Cucuta", "Bucaramanga", "Pereira", "Santa Marta", "Ibague", 
                      "Pasto", "Manizales", "Neiva")
n_ciudades <- length(ciudades_nombres)

Agregamos datos necesarios para graficar la ruta:

coords <- data.frame(
  ciudad = ciudades_nombres,
  lon = c(-74.0721, -76.5225, -75.5736, -74.7964, -75.5144, -72.5074, -73.1198, -75.6961, -74.2050, -75.2174, -77.2811, -75.4849, -75.2819),
  lat = c(4.7110, 3.4372, 6.2442, 10.9639, 10.3910, 7.8939, 7.1139, 4.8133, 11.2404, 4.4389, 1.2136, 5.0689, 2.9273)
)

colombia <- ne_countries(scale = "medium", country = "Colombia", returnclass = "sf")

Definimos parámetros de costos iniciales para el salario, consumo de combustible y precio de gasolina:

salario_hora <- 6189  # COP/hora
consumo_combustible <- 60  # km/galón
precio_gasolina <- 15827  # COP/galón

Calculamos el costo de la gasolina, el tiempo de los viajes y el costo del salario:

df_costo_gasolina <- (df_distancias/consumo_combustible) * precio_gasolina
df_costo_salario <- df_tiempos_viajes * salario_hora

Ahora calculamos el costo total de desplazamiento entre las ciudades:

df_costo_total <- df_costo_gasolina + df_costo_salario + df_peajes
df_costo_total

##                Bogota     Cali Medellin Barranquilla Cartagena   Cucuta
## Bogota            0.0 361925.4 262810.7    584676.17 614344.66 301176.8
## Cali         356731.7      0.0 330387.9    779820.96 783989.45 568161.5
## Medellin     262908.4 330379.0      0.0    461341.95 456008.71 348705.3
## Barranquilla 610423.0 785956.8 466623.7         0.00  93896.97 406218.1
## Cartagena    630503.6 786109.9 456054.5     88974.36      0.00 432498.1
## Cucuta       299859.5 583224.4 346691.3    403964.37 435006.82      0.0
## Bucaramanga  242248.0 458501.1 221966.2    377287.57 408330.02 124826.5
## Pereira      259752.8 145069.2 184248.8    664381.57 622450.05 511174.9
## Santa Marta  533350.8 697966.9 461433.8     84044.68 182559.96 329145.9
## Ibague       168694.8 181424.3 267546.1    572620.39 602288.88 391939.6
## Pasto        529225.7 192144.7 502781.9    961714.95 959682.17 830155.5
## Manizales    146148.3 196551.3 160953.3    641086.42 617754.91 381661.9
## Neiva        228458.7 199599.9 355108.2    676971.93 706640.42 496291.2
##              Bucaramanga   Pereira Santa Marta    Ibague    Pasto Manizales
## Bogota          244064.3 264154.77   523767.01 167357.68 520448.5 146073.68
## Cali            441421.7 144583.40   696403.66 175742.21 193372.9 197047.14
## Medellin        221966.9 184095.29   460724.65 261494.23 504104.8 160821.96
## Barranquilla    381176.5 670374.81    85071.98 578606.93 969182.6 647101.49
## Cartagena       407456.2 624426.24   184084.41 603037.79 963134.7 619752.65
## Cucuta          123083.6 526542.41   321638.12 409450.99 846550.2 409967.37
## Bucaramanga          0.0 404517.39   294961.32 287425.96 724526.9 287944.06
## Pereira         387136.6      0.00   642565.99  78663.35 318795.0  53822.79
## Santa Marta     304104.0 643983.19        0.00 561717.03 963992.7 527409.87
## Ibague          267901.6  83553.60   511711.49      0.00 371448.3 153420.41
## Pasto           706115.4 316977.39   961097.65 364436.20      0.0 369441.13
## Manizales       257623.6  53819.93   525969.12 147518.87 370277.1      0.00
## Neiva           372253.1 206229.22   616064.75 125731.86 233667.5 259796.03
##                 Neiva
## Bogota       228171.4
## Cali         193884.5
## Medellin     354058.6
## Barranquilla 682458.7
## Cartagena    706887.9
## Cucuta       508914.9
## Bucaramanga  391277.7
## Pereira      201305.9
## Santa Marta  665567.1
## Ibague       126547.7
## Pasto        231941.9
## Manizales    253861.2
## Neiva             0.0

#exportar dtf_costo_total a csv, verificando que no exista el archivo
if (file.exists("df_costo_total.csv")) {
  file.remove("df_costo_total.csv")
  write.csv(df_costo_total, "df_costo_total.csv")
} else {
  write.csv(df_costo_total, "df_costo_total.csv")
}

Ahora, dado este dataset de costos, utilizaremos tanto el algoritmo de colonia de hormigas y un algoritmo evolutivo para encontrar el orden correcto.

Algoritmo de Optimización por Colonia de Hormigas (ACO)

La Optimización por Colonia de Hormigas (ACO) es “una técnica probabilística para solucionar problemas computacionales que pueden reducirse a buscar los mejores caminos o rutas en grafos”, como se describe en la literatura (Wikipedia, s.f.). Su funcionamiento se inspira directamente en cómo las hormigas encuentran eficientemente fuentes de alimento: depositan feromonas en sus trayectos, y las rutas más efectivas acumulan mayores concentraciones de esta sustancia, atrayendo a más hormigas.

El algoritmo ACO simula este comportamiento mediante:

Hormigas Artificiales: Agentes que construyen soluciones candidatas (rutas).

Feromonas Virtuales: Información asociada a los componentes de la solución (e.g., las conexiones entre ciudades) que se actualiza según la calidad de las soluciones encontradas. Rutas de mejor calidad refuerzan más intensamente sus componentes.

Selección Probabilística: Las hormigas eligen los siguientes pasos basándose tanto en la intensidad de la feromona como, opcionalmente, en información heurística (e.g., la distancia directa), favoreciendo caminos prometedores pero permitiendo la exploración.

A continuación, se presenta la función ant_colony_optimization_r_con_historial desarrollada en R para implementar el algoritmo de optimización por colonia de hormigas. Esta función está diseñada para resolver el Problema del Viajante de Comercio (TSP), buscando la ruta de menor costo total. Además de encontrar la mejor ruta, registra el historial de la mejor solución encontrada en cada iteración, lo cual es útil para analizar la convergencia y visualizar el proceso.

# --- Función ACO Única (con historial detallado de rutas) ---
# Nota: Esta función devuelve tanto el resultado final como el historial.
ant_colony_optimization_r_con_historial <- function(costos, num_ants = 10, num_iterations = 100, alpha = 1.0, beta = 2.0, rho = 0.5, Q = 100, verbose = TRUE) {

  # Convertir a matriz y validar
  cost_matrix <- as.matrix(costos)
  n_nodes <- ncol(cost_matrix)
  if (nrow(cost_matrix) != n_nodes) stop("La matriz de costos debe ser cuadrada.")

  # Calcular visibilidad (eta)
  visibility_matrix <- cost_matrix
  diag(visibility_matrix) <- Inf # No viajar de un nodo a sí mismo
  eta <- 1 / visibility_matrix
  eta[is.infinite(eta)] <- 0 # Manejar división por cero o Inf si el costo era 0 o Inf
  eta[visibility_matrix == 0] <- 0 # Asegurar que costo 0 no de visibilidad infinita

  # Inicializar feromona
  mean_cost <- mean(cost_matrix[cost_matrix > 0 & !is.infinite(cost_matrix)], na.rm = TRUE)
  initial_pheromone <- if (is.na(mean_cost) || mean_cost <= 0 || !is.finite(mean_cost)) 1.0 else 1 / (n_nodes * mean_cost)
  pheromone_matrix <- matrix(initial_pheromone, nrow = n_nodes, ncol = n_nodes)
  diag(pheromone_matrix) <- 0 # Sin feromona en el propio nodo

  # Variables para guardar la mejor solución encontrada
  best_tour <- NULL
  best_length <- Inf

  # Estructuras para el historial
  history_best_length <- numeric(num_iterations) # Guarda el mejor costo hasta esa iteración
  history_best_tours <- vector("list", num_iterations) # Guarda la mejor ruta hasta esa iteración

  # --- Bucle principal de iteraciones ---
  for (iter in 1:num_iterations) {
    all_tours <- vector("list", num_ants) # Rutas de todas las hormigas en esta iteración
    all_lengths <- numeric(num_ants)      # Costos de todas las rutas en esta iteración
    iter_found_new_best <- FALSE           # Flag para indicar si se mejoró en esta iteración

    # --- Bucle de hormigas ---
    for (ant in 1:num_ants) {
      current_tour <- numeric(n_nodes) # Ruta de la hormiga actual
      visited <- logical(n_nodes)      # Nodos visitados por la hormiga actual

      # Empezar en un nodo aleatorio
      start_node <- sample(1:n_nodes, 1)
      current_tour[1] <- start_node
      visited[start_node] <- TRUE
      current_node <- start_node

      # --- Construcción de la ruta de la hormiga ---
      for (step in 2:n_nodes) {
        possible_nodes <- which(!visited) # Nodos aún no visitados

        if (length(possible_nodes) == 0) break # Si no hay más nodos, terminar (raro en TSP)

        # Obtener feromonas y visibilidad hacia los nodos posibles
        pheromones_to_possible <- pheromone_matrix[current_node, possible_nodes]
        eta_to_possible <- eta[current_node, possible_nodes]

        # Manejar valores no finitos o muy pequeños para evitar NaN/Inf en probabilidades
        pheromones_to_possible[!is.finite(pheromones_to_possible)] <- 1e-10
        eta_to_possible[!is.finite(eta_to_possible)] <- 1e-10
        pheromones_to_possible[pheromones_to_possible < 1e-10] <- 1e-10
        eta_to_possible[eta_to_possible < 1e-10] <- 1e-10

        # Calcular numerador de la probabilidad de transición
        prob_numerator <- (pheromones_to_possible^alpha) * (eta_to_possible^beta)
        sum_prob <- sum(prob_numerator)

        # Calcular probabilidades finales
        probabilities <- if (sum_prob > 1e-10 && is.finite(sum_prob)) {
          prob_numerator / sum_prob
        } else {
          # Si suma es cero o no finita, asignar probabilidad uniforme
          rep(1 / length(possible_nodes), length(possible_nodes))
        }

        # Asegurar que las probabilidades sumen 1 (manejo de errores numéricos)
        if(any(!is.finite(probabilities)) || sum(probabilities) < 1e-10) {
           probabilities <- rep(1 / length(possible_nodes), length(possible_nodes))
        }
        probabilities <- probabilities / sum(probabilities) # Normalizar

        # Elegir el siguiente nodo basado en las probabilidades
        chosen_index <- if (length(possible_nodes) == 1) {
          1 # Si solo queda uno, elegir ese
        } else {
          sample(1:length(possible_nodes), 1, prob = probabilities)
        }
        next_node <- possible_nodes[chosen_index]

        # Actualizar ruta y estado de visitados
        current_tour[step] <- next_node
        visited[next_node] <- TRUE
        current_node <- next_node
      } # Fin construcción de ruta de la hormiga (step)

      # --- Calcular costo de la ruta construida ---
      current_length <- 0
      valid_tour <- TRUE
      # Verificar si se visitaron todos los nodos (importante si hubo 'break')
      if (!all(1:n_nodes %in% current_tour)) {
        valid_tour <- FALSE
        current_length <- Inf
      } else {
        # Calcular costo sumando segmentos
        for (i in 1:(n_nodes - 1)) {
          from_node <- current_tour[i]
          to_node <- current_tour[i+1]
          cost_segment <- cost_matrix[from_node, to_node]
          if (!is.finite(cost_segment)) { # Ruta inválida si un segmento tiene costo Inf
            current_length <- Inf
            break
          }
          current_length <- current_length + cost_segment
        }
        # Añadir costo de regreso al nodo inicial si la ruta es finita hasta ahora
        if (is.finite(current_length)) {
          from_node <- current_tour[n_nodes]
          to_node <- current_tour[1]
          cost_segment <- cost_matrix[from_node, to_node]
          if (!is.finite(cost_segment)) {
            current_length <- Inf # Ruta inválida si el regreso es Inf
          } else {
            current_length <- current_length + cost_segment
          }
        }
      } # Fin cálculo de costo

      # Guardar ruta y costo de la hormiga
      all_tours[[ant]] <- if(valid_tour && is.finite(current_length)) current_tour else NULL
      all_lengths[ant] <- current_length

      # Actualizar la mejor solución global si esta hormiga encontró una mejor
      if (is.finite(current_length) && current_length < best_length) {
        if(verbose) cat(sprintf("Iter %d, Hormiga %d: Nueva mejor longitud = %.2f\n", iter, ant, current_length))
        best_length <- current_length
        best_tour <- current_tour
        iter_found_new_best <- TRUE
      }
    } # Fin bucle de hormigas (ant)

    # --- Actualización de Feromonas ---
    # 1. Evaporación
    pheromone_matrix <- (1 - rho) * pheromone_matrix

    # 2. Deposición
    delta_pheromone <- matrix(0, nrow = n_nodes, ncol = n_nodes) # Cambios en feromona para esta iteración
    for (ant in 1:num_ants) {
      ant_tour <- all_tours[[ant]]
      ant_length <- all_lengths[ant]
      # Depositar solo si la ruta es válida, completa y finita
      if (!is.null(ant_tour) && length(ant_tour) == n_nodes && all(ant_tour > 0) && is.finite(ant_length) && ant_length > 0) {
        deposit_amount <- Q / ant_length # Cantidad de feromona a depositar
        for (i in 1:(n_nodes - 1)) {
          node1 <- ant_tour[i]
          node2 <- ant_tour[i+1]
          # Depositar en ambas direcciones (simétrico para TSP)
          delta_pheromone[node1, node2] <- delta_pheromone[node1, node2] + deposit_amount
          delta_pheromone[node2, node1] <- delta_pheromone[node2, node1] + deposit_amount
        }
        # Depositar en el segmento de cierre (último a primero)
        node1 <- ant_tour[n_nodes]
        node2 <- ant_tour[1]
        delta_pheromone[node1, node2] <- delta_pheromone[node1, node2] + deposit_amount
        delta_pheromone[node2, node1] <- delta_pheromone[node2, node1] + deposit_amount
      }
    }
    pheromone_matrix <- pheromone_matrix + delta_pheromone

    # Asegurar un mínimo de feromona y cero en diagonal
    pheromone_matrix[pheromone_matrix < 1e-10] <- 1e-10
    diag(pheromone_matrix) <- 0

    # --- Guardar historial de la iteración ---
    # Guarda el mejor costo encontrado *hasta ahora* en esta iteración
    history_best_length[iter] <- best_length
    # Guarda la mejor ruta encontrada *hasta ahora* en esta iteración
    history_best_tours[[iter]] <- best_tour

    # Mensaje de progreso periódico
    if(verbose && iter %% 10 == 0) {
      cat(sprintf("--- Fin Iteración %d: Mejor costo hasta ahora = %.2f ---\n", iter, best_length))
    }

  } # Fin bucle principal de iteraciones (iter)

  # Mensaje final
  if (verbose && is.infinite(best_length)) {
    cat("ACO: No se encontró ninguna ruta válida finita.\n")
  } else if (verbose) {
    cat(sprintf("ACO: Optimización finalizada. Mejor longitud encontrada = %.2f\n", best_length))
  }

  # --- Preparar y devolver resultados ---
  # Intentar obtener nombres de ciudades si 'ciudades_nombres' existe y best_tour no es NULL
  best_tour_cities_result <- NULL
  if (!is.null(best_tour) && exists("ciudades_nombres")) {
     best_tour_cities_result <- try(ciudades_nombres[best_tour], silent = TRUE)
     if (inherits(best_tour_cities_result, "try-error")) {
         best_tour_cities_result <- "Error al obtener nombres" # O asignar NULL
     }
  } else if (is.null(best_tour)){
      best_tour_cities_result <- "No disponible (sin ruta)"
  } else {
      best_tour_cities_result <- "No disponible (ciudades_nombres no existe)"
  }


  return(list(
    best_tour_indices = best_tour,          # Vector de índices de la mejor ruta
    best_tour_cost = best_length,           # Costo numérico de la mejor ruta
    best_tour_cities = best_tour_cities_result, # Vector de nombres de ciudades (o mensaje)
    convergence_cost = history_best_length, # Vector numérico de costos por iteración
    history_best_tours = history_best_tours # Lista de vectores de índices (rutas) por iteración
  ))
}

Ahora, ejecutamos la función ant_colony_optimization_r_con_historial con un conjunto de parámetros definidos. Se establece una semilla (set.seed) para garantizar la reproducibilidad de los resultados, dado el componente probabilístico del algoritmo.

Adicional, en la figura 23 se observa como el algoritmo va obtendiendo mejores resultados al optimizar la función.

# Parámetros de ejecución
num_ants_aco <- 30
num_iterations_aco <- 100

set.seed(456) # Semilla para reproducibilidad
resultado_aco <- ant_colony_optimization_r_con_historial(
  costos = df_costo_total, # Asume que df_costo_total existe
  num_ants = num_ants_aco,
  num_iterations = num_iterations_aco,
  alpha = 1.0,
  beta = 3.0,
  rho = 0.3,
  Q = 100,
  verbose = TRUE
)

## Iter 1, Hormiga 1: Nueva mejor longitud = 6242355.38
## Iter 1, Hormiga 3: Nueva mejor longitud = 6129624.32
## Iter 1, Hormiga 4: Nueva mejor longitud = 3691924.90
## Iter 9, Hormiga 14: Nueva mejor longitud = 3521789.45
## --- Fin Iteración 10: Mejor costo hasta ahora = 3521789.45 ---
## --- Fin Iteración 20: Mejor costo hasta ahora = 3521789.45 ---
## Iter 26, Hormiga 12: Nueva mejor longitud = 3305521.55
## Iter 28, Hormiga 22: Nueva mejor longitud = 3299687.58
## --- Fin Iteración 30: Mejor costo hasta ahora = 3299687.58 ---
## --- Fin Iteración 40: Mejor costo hasta ahora = 3299687.58 ---
## --- Fin Iteración 50: Mejor costo hasta ahora = 3299687.58 ---
## --- Fin Iteración 60: Mejor costo hasta ahora = 3299687.58 ---
## --- Fin Iteración 70: Mejor costo hasta ahora = 3299687.58 ---
## Iter 80, Hormiga 14: Nueva mejor longitud = 2860689.07
## --- Fin Iteración 80: Mejor costo hasta ahora = 2860689.07 ---
## --- Fin Iteración 90: Mejor costo hasta ahora = 2860689.07 ---
## --- Fin Iteración 100: Mejor costo hasta ahora = 2860689.07 ---
## ACO: Optimización finalizada. Mejor longitud encontrada = 2860689.07

print("--- Resultado Final - Ejecución ACO ---")

## [1] "--- Resultado Final - Ejecución ACO ---"

if (!is.null(resultado_aco) && is.finite(resultado_aco$best_tour_cost)) {
    print(paste("Mejor costo encontrado:", format(resultado_aco$best_tour_cost, scientific = FALSE, big.mark=",")))
    # Verificar si best_tour_cities existe y no es NULL antes de imprimir
    # Asume que 'ciudades_nombres' fue usado correctamente dentro de la función ACO o existe globalmente
    if (!is.null(resultado_aco$best_tour_cities) && !inherits(resultado_aco$best_tour_cities, "try-error") && !startsWith(resultado_aco$best_tour_cities[1], "Error") && !startsWith(resultado_aco$best_tour_cities[1], "No disponible")) {
       print(paste("Mejor ruta (ciudades):", paste(resultado_aco$best_tour_cities, collapse = " -> ")))
    } else if (!is.null(resultado_aco$best_tour_indices)) {
        print(paste("Mejor ruta (índices):", paste(resultado_aco$best_tour_indices, collapse=" -> ")))
        print("Mejor ruta (ciudades): No disponible o error al obtener nombres.")
    } else {
       print("Mejor ruta (ciudades e índices): No disponible (best_tour fue NULL o inválido).")
    }

    # --- Modificación Aquí: Graficar convergencia del costo CON CAPTION ---
    tryCatch({
        # --- Inicio: Ajustar y guardar márgenes ---
        old_mar <- par("mar") # Guardar márgenes originales [bottom, left, top, right]
        # Incrementar margen inferior (el primer valor) para dar espacio al caption
        new_mar <- old_mar + c(2, 0, 0, 0) # Añadir 2 líneas al margen inferior
        par(mar = new_mar)
        # --- Fin: Ajustar y guardar márgenes ---

        # 1. Generar el gráfico con títulos y etiquetas
        plot(1:num_iterations_aco, resultado_aco$convergence_cost, type='l',
             col="steelblue", lwd=2, # Un poco de estilo
             xlab="Iteración",
             ylab="Mejor Costo Encontrado",
             main="Convergencia del Costo en ACO")
        grid() # Añadir una rejilla

        # 2. Añadir el caption en el margen inferior
        caption_text <- "Figura 23: Evolución del mejor costo acumulado encontrado por ACO en cada iteración."
        # Ajusta 'line' (e.g., 4, 4.5, 5) si el caption se superpone con xlab
        mtext(caption_text, side = 1, line = 5, adj = 0, cex = 0.8)

        # --- Inicio: Restaurar márgenes ---
        par(mar = old_mar) # Restaurar márgenes originales para no afectar otros gráficos
        # --- Fin: Restaurar márgenes ---

      },
      error=function(e) {
        print(paste("Error al graficar convergencia ACO:", e$message))
        # Asegurarse de restaurar márgenes incluso si hay error y old_mar existe
        if(exists("old_mar") && !is.null(old_mar)) {
           try(par(mar = old_mar), silent = TRUE) # Intenta restaurar
        }
      }
    )
    # --- Fin Modificación ---

} else {
    print("ACO no encontró una solución válida o la ejecución falló.")
    # Podrías imprimir resultado_aco aquí para depurar si es necesario
    # print(resultado_aco)
}

## [1] "Mejor costo encontrado: 2,860,689"
## [1] "Mejor ruta (ciudades): Pasto -> Neiva -> Ibague -> Manizales -> Pereira -> Bogota -> Bucaramanga -> Cucuta -> Barranquilla -> Cartagena -> Santa Marta -> Medellin -> Cali"

Para entender mejor cómo el algoritmo refina la solución, generamos una animación que muestra la mejor ruta encontrada hasta cada iteración. Primero, preparamos los datos del historial de rutas, convirtiendo los índices de las rutas en coordenadas geográficas y estructurándolos para gganimate .

# --- Preparar Datos del Historial para gganimate ---
# Asume que 'resultado_aco' (con historial), 'coords', y 'n_ciudades' ya existen

historial_rutas <- resultado_aco$history_best_tours
historial_costos <- resultado_aco$convergence_cost
all_segments_history <- data.frame() # DataFrame para segmentos históricos

# Validar que el historial existe y tiene elementos
if (!is.null(historial_rutas) && length(historial_rutas) > 0) {
    
    # Usar el número real de iteraciones guardadas
    num_iter_guardadas <- length(historial_rutas) 
    
    for (iter in 1:num_iter_guardadas) { 
      current_tour_indices <- historial_rutas[[iter]]
      # Asegurarse de que el índice de costo no exceda la longitud real
      current_cost <- if(iter <= length(historial_costos)) historial_costos[iter] else NA 
      
      # Procesar solo si hay una ruta válida y costo válido
      if (!is.null(current_tour_indices) && length(current_tour_indices) == n_ciudades && all(current_tour_indices > 0) && is.finite(current_cost)) {
        
        ruta_iter_coords <- coords[current_tour_indices, ]
        ruta_iter_coords <- rbind(ruta_iter_coords, ruta_iter_coords[1, ]) # Cerrar ciclo
        
        segmentos_iter <- data.frame(
          lon_inicio = head(ruta_iter_coords$lon, -1),
          lat_inicio = head(ruta_iter_coords$lat, -1),
          lon_fin = tail(ruta_iter_coords$lon, -1),
          lat_fin = tail(ruta_iter_coords$lat, -1),
          iteracion = iter, 
          costo = current_cost 
        )
        all_segments_history <- rbind(all_segments_history, segmentos_iter)
      } 
    }
}

# Verificar si se generaron segmentos
if(nrow(all_segments_history) == 0) {
    print("ADVERTENCIA: No se pudieron generar datos de segmentos del historial de ACO para la animación.")
} else {
    print(paste("Se procesaron datos de", length(unique(all_segments_history$iteracion)), "iteraciones de ACO para la animación."))
}

## [1] "Se procesaron datos de 100 iteraciones de ACO para la animación."

# --- Animación del Proceso de Mejora de ACO ---
# Asume que 'all_segments_history', 'colombia', y 'coords' existen y son válidos.
# También asume que las librerías ggplot2, gganimate, sf, gifski están cargadas.

# Solo intentar animar si hay datos válidos
if(exists("all_segments_history") && nrow(all_segments_history) > 0) {

    # Crear mapa base
    mapa_base_historia <- ggplot() +
      geom_sf(data = colombia, fill = "ivory", color = "gray50") +
      geom_point(data = coords, aes(x = lon, y = lat), color = "darkblue", size = 3.0, alpha = 0.6) +
      theme_minimal(base_size = 11) +
      theme(plot.title = element_text(hjust = 0.5, face="bold"),
            plot.subtitle = element_text(hjust = 0.5)) +
      coord_sf(xlim = range(coords$lon) + c(-1.2, 1.2), ylim = range(coords$lat) + c(-1.2, 1.2), expand = FALSE)
    
    # Crear objeto de animación
    animacion_proceso <- mapa_base_historia +
      geom_segment(data = all_segments_history,
                   aes(x = lon_inicio, y = lat_inicio, xend = lon_fin, yend = lat_fin),
                   color = "darkorange", size = 1.0, alpha = 0.8) + 
      transition_manual(iteracion) +  

      labs(
        title = "Proceso de Optimización ACO - Iteración {frame}", 
        subtitle = "Mostrando la mejor ruta encontrada hasta esta iteración", 
        caption = "Figura 24: Proceso de optimización ACO - Iteraciones. (Fuente: Elaboración propia, 2025).", y = "Latitud"
      ) +
      ease_aes('linear') 
    
    # Renderizar animación
    fps_proceso <- 5 # Velocidad de la animación
    # Usar el número de frames únicos disponibles en los datos procesados
    total_frames_proceso <- length(unique(all_segments_history$iteracion)) 
    
    print(paste("Generando animación del proceso ACO con", total_frames_proceso, "frames a", fps_proceso, "fps."))
    
    # Usar tryCatch para manejar posibles errores durante la renderización
    anim_render_proceso <- tryCatch({
        animate(
          animacion_proceso, 
          nframes = total_frames_proceso, 
          fps = fps_proceso, 
          width = 800, 
          height = 700, 
          renderer = gifski_renderer()
        )
    }, error = function(e) {
        print(paste("Error durante la renderización de la animación:", e$message))
        return(NULL) # Devuelve NULL si falla
    })
    
    # Solo proceder si la renderización fue exitosa
    if (!is.null(anim_render_proceso)) {
      # Visualizar en Rmd/Viewer (se imprime automáticamente por defecto en Rmd si es la última expresión)
      # print(anim_render_proceso) # Descomentar si no se muestra automáticamente
      
      # Guardar
      tryCatch({
          anim_save("proceso_optimizacion_aco_optimizado.gif", animation = anim_render_proceso)
          print("Animación del proceso ACO guardada como proceso_optimizacion_aco_optimizado.gif")
      }, error = function(e) {
          print(paste("Error al guardar la animación:", e$message))
      })
      
    } # Fin if !is.null(anim_render_proceso)

} else {
  print("No hay datos históricos válidos para generar la animación del proceso ACO.")
}

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## [1] "Generando animación del proceso ACO con 100 frames a 5 fps."

## Warning in device(files[i], width = 800, height = 700, units = "in", res =
## 192): cairo error 'invalid value (typically too big) for the size of the input
## (surface, pattern, etc.)'

## [1] "Error durante la renderización de la animación: unable to start device 'png'"

La animación resultante (Figura 24) ilustra cómo la mejor ruta encontrada por el algoritmo ACO evoluciona a lo largo de las iteraciones, tendiendo a estabilizarse a medida que el algoritmo converge hacia una solución de bajo costo.

knitr::include_graphics("proceso_optimizacion_aco_optimizado.gif")

Continuando con con el Algoritmo Genético:

Continuando con las metaheurísticas, ahora abordamos los Algoritmos Genéticos (AG), una clase de algoritmos de optimización inspirados en el proceso de la evolución natural y la genética. Los AG son particularmente efectivos para problemas complejos donde el espacio de búsqueda es grande, como el Problema del Viajante de Comercio.

El enfoque se basa en mantener y evolucionar una población de soluciones candidatas al problema. En nuestro contexto del TSP, cada ‘individuo’ de la población es una ruta o secuencia específica de visita a las ciudades. Cada ruta se representa mediante una codificación, análoga a un cromosoma (usualmente, una permutación de los índices o nombres de las ciudades).

El algoritmo opera en ciclos llamados generaciones. En cada generación, se simula el proceso evolutivo mediante los siguientes pasos clave [inspirados en conceptos descritos por Wikipedia (s.f.) y otros]:

Evaluación: Se calcula la aptitud (fitness) de cada ruta en la población. Para el TSP, una mayor aptitud corresponde a un menor costo total del recorrido (usualmente se define como el inverso del costo o alguna función decreciente del costo).
Selección: Se seleccionan las rutas ‘padres’ que contribuirán a la siguiente generación. Las rutas con mayor aptitud tienen una mayor probabilidad de ser elegidas, simulando la “supervivencia del más apto”.
Cruzamiento (Crossover): Se combinan pares de rutas padres para crear nuevas rutas ‘hijas’. Este proceso mezcla segmentos o características de las rutas progenitoras, buscando generar nuevas combinaciones potencialmente mejores.
Mutación: Se aplican pequeñas alteraciones aleatorias a algunas de las nuevas rutas (o a veces a individuos existentes). En el TSP, esto podría implicar intercambiar la posición de dos ciudades en la secuencia. La mutación introduce diversidad genética en la población y ayuda a evitar que el algoritmo se estanque prematuramente en soluciones subóptimas.
Reemplazo: Se forma la población para la siguiente generación, decidiendo qué individuos (padres, hijos) sobreviven. Existen diversas estrategias, como reemplazar a los menos aptos o mantener a los mejores (‘elitismo’).

A través de la aplicación repetida de estos operadores genéticos (selección, cruzamiento, mutación) durante muchas generaciones, la población tiende a converger hacia rutas de muy alta calidad (bajo costo).

A continuación, se presenta la función en R que implementa el Algoritmo Genético diseñado para resolver nuestro problema del viajante.

calculate_tour_cost <- function(tour, cost_matrix) {
  total_cost <- 0
  n <- length(tour)
  
  if (n < 2) return(Inf) 
  
  for (i in 1:(n - 1)) { 
    idx1 <- tour[i]
    idx2 <- tour[i+1]
    
    # Validación de índices
    if(any(is.na(c(idx1,idx2))) || idx1 < 1 || idx1 > nrow(cost_matrix) || idx2 < 1 || idx2 > ncol(cost_matrix)) {
      return(Inf)
    }
      
    cost <- cost_matrix[idx1, idx2]
    if (!is.finite(cost)) {
      return(Inf)
    }
      
    total_cost <- total_cost + cost 
  }
  
  
  idx1 <- tour[n]
  idx2 <- tour[1]
  
  if(any(is.na(c(idx1,idx2))) || idx1 < 1 || idx1 > nrow(cost_matrix) || idx2 < 1 || idx2 > ncol(cost_matrix)) {
    return(Inf)
  }
    
  cost_return <- cost_matrix[idx1, idx2]
  if (!is.finite(cost_return)) {
    return(Inf)
  }
    
  total_cost <- total_cost + cost_return
  return(total_cost)
}

fitness_function <- function(tour, cost_matrix) { 
  cost <- calculate_tour_cost(tour, cost_matrix)
  if (is.infinite(cost)) { 
    # Devolver valor muy negativo pero finito si el costo es Infinito
    return(-1e+100) 
  } else { 
    # Devolver negativo del costo
    return(-cost) 
  } 
}

# Listas externas para guardar historial (se deben limpiar antes de llamar a ga)
ga_history_best_tours <- list()
ga_history_best_costs <- list()

# Función Monitor personalizada para GA
monitor_ga_history <- function(object, ...) { 
  iter <- object@iter
  best_fitness_current_gen <- max(object@fitness, na.rm = TRUE)
  
  # No hacer nada si no hay fitness finito
  if (!is.finite(best_fitness_current_gen)) {
    return() 
  }
    
  best_indices_current_gen <- which(object@fitness == best_fitness_current_gen)
  
  if (length(best_indices_current_gen) > 0) { 
    # Tomar la primera mejor solución en caso de empate
    best_solution_current_gen <- object@population[best_indices_current_gen[1], ]
    actual_cost <- -best_fitness_current_gen
    
    # Guardar en listas externas usando superasignación (<<-)
    ga_history_best_tours[[as.character(iter)]] <<- as.vector(best_solution_current_gen)
    ga_history_best_costs[[as.character(iter)]] <<- actual_cost
    
    # Imprimir progreso cada 20 generaciones o en la primera
    if (iter %% 20 == 0 || iter == 1) {
      cat(sprintf("GA Gen %d: Mejor Costo (Gen actual) = %.2f\n", iter, actual_cost)) 
    }
  } # Fin if best_indices
} # Fin monitor_ga_history



# Resetear listas de historial antes de la ejecución
ga_history_best_tours <- list()
ga_history_best_costs <- list() 

# Fijar semilla y número de generaciones
set.seed(456)
num_generations_ga <- 200 

# Ejecutar GA usando tryCatch para manejo de errores
resultado_ga <- tryCatch({
  ga(
    type = "permutation", 
    fitness = fitness_function, 
    cost_matrix = df_costo_total, 
    lower = 1, 
    upper = n_ciudades, 
    popSize = 100, 
    maxiter = num_generations_ga, 
    run = 100, 
    pmutation = 0.2, 
    pcrossover = 0.8, 
    elitism = max(1, floor(0.1 * 100)), 
    monitor = monitor_ga_history, # Usar el monitor personalizado
    seed = 123
  )
}, error = function(e) {
    # Imprimir error y devolver NULL si falla
    print(paste("Error durante la ejecución de GA:", e$message))
    return(NULL)
})

## GA Gen 1: Mejor Costo (Gen actual) = 3489938.86
## GA Gen 20: Mejor Costo (Gen actual) = 3007253.95
## GA Gen 40: Mejor Costo (Gen actual) = 2661123.84
## GA Gen 60: Mejor Costo (Gen actual) = 2661123.84
## GA Gen 80: Mejor Costo (Gen actual) = 2453028.09
## GA Gen 100: Mejor Costo (Gen actual) = 2453028.09
## GA Gen 120: Mejor Costo (Gen actual) = 2452588.14
## GA Gen 140: Mejor Costo (Gen actual) = 2452588.14
## GA Gen 160: Mejor Costo (Gen actual) = 2452588.14
## GA Gen 180: Mejor Costo (Gen actual) = 2403977.63
## GA Gen 200: Mejor Costo (Gen actual) = 2403977.63

# Verificar resultado de la ejecución
if (!is.null(resultado_ga)) { 
  print("Ejecución GA finalizada.")
  print(paste("Historial GA guardado para", length(ga_history_best_tours), "generaciones.")) 
} else { 
  print("Ejecución GA falló.")
}

## [1] "Ejecución GA finalizada."
## [1] "Historial GA guardado para 200 generaciones."

print("--- Resultados Finales GA ---")

## [1] "--- Resultados Finales GA ---"

# Solo mostrar si la ejecución fue exitosa
if (!is.null(resultado_ga)) {
  
  # Inicializar variables de resultado
  best_tour_indices_ga <- NULL
  best_tour_cost_ga <- Inf
  best_tour_cities_ga <- "No disponible"
  
  # Extraer la mejor solución final si es válida
  if (length(resultado_ga@solution) > 0 && is.finite(resultado_ga@fitnessValue)) { 
    best_tour_indices_ga <- resultado_ga@solution[1,] # Primera mejor solución
    best_tour_cost_ga <- -resultado_ga@fitnessValue  # Costo real (positivo)
    
    # Obtener nombres de ciudades (si existen)
    best_tour_cities_ga <- tryCatch(
        ciudades_nombres[best_tour_indices_ga], 
        error = function(e) "Error al obtener nombres"
    ) 
  } else { 
    print("Advertencia: No se encontró solución final GA válida.")
  }
  
  # Imprimir resultados finales formateados
  print(paste("Mejor costo GA:", format(best_tour_cost_ga, scientific = FALSE, big.mark=",")))
  if (!is.null(best_tour_indices_ga)) {
     print(paste("Mejor ruta GA (índices):", paste(best_tour_indices_ga, collapse = " -> "))) 
  } else {
      print("Mejor ruta GA (índices): No disponible")
  }
  print(paste("Mejor ruta GA (ciudades):", paste(best_tour_cities_ga, collapse = " -> ")))
  
  # Graficar convergencia estándar del paquete GA
 tryCatch({
      # --- Inicio: Ajustar y guardar márgenes ---
      old_mar <- par("mar") # Guardar márgenes originales [bottom, left, top, right]
      # Incrementar margen inferior (el primer valor) para dar espacio al caption
      # Aumentamos de 5.1 (default usual) a 6.1 o 7.1, probemos con +2 líneas
      new_mar <- old_mar + c(2, 0, 0, 0)
      par(mar = new_mar)
      # --- Fin: Ajustar y guardar márgenes ---

      # 1. Generar el gráfico base
      plot(resultado_ga)

      # 2. Añadir las etiquetas y el título
      title(main = "Convergencia del Algoritmo Genético (GA)")

      caption_text <- "Figura 25: Evolución del mejor fitness (verde) y fitness promedio (azul) por generación."

      mtext(caption_text, side = 1, line = 5, adj = 0, cex = 0.8)

      # --- Inicio: Restaurar márgenes ---
      par(mar = old_mar) # Restaurar márgenes originales para no afectar otros gráficos
      # --- Fin: Restaurar márgenes ---

    },
    error=function(e) {
      print(paste("Error al graficar o rotular GA:", e$message))
      # Asegurarse de restaurar márgenes incluso si hay error y old_mar existe
      if(exists("old_mar") && !is.null(old_mar)) {
         try(par(mar = old_mar), silent = TRUE) # Intenta restaurar
      }
    }
  )
  
} else { 
  print("No hay resultados finales GA para mostrar.") 
}

## [1] "Mejor costo GA: 2,403,978"
## [1] "Mejor ruta GA (índices): 3 -> 12 -> 8 -> 2 -> 11 -> 13 -> 10 -> 1 -> 7 -> 6 -> 9 -> 4 -> 5"
## [1] "Mejor ruta GA (ciudades): Medellin -> Manizales -> Pereira -> Cali -> Pasto -> Neiva -> Ibague -> Bogota -> Bucaramanga -> Cucuta -> Santa Marta -> Barranquilla -> Cartagena"

# Data frame para almacenar todos los segmentos de todas las generaciones
all_segments_history_ga <- data.frame()

# Procesar solo si hay historial y coordenadas
if (length(ga_history_best_tours) > 0 && exists("coords")) {
  
  generaciones_guardadas <- names(ga_history_best_tours)
  
  # Iterar sobre cada generación guardada
  for (gen_char in generaciones_guardadas) {
    gen <- as.numeric(gen_char)
    current_tour_indices <- ga_history_best_tours[[gen_char]]
    current_cost <- ga_history_best_costs[[gen_char]]
    
    # Validar ruta y costo antes de procesar
    if (!is.null(current_tour_indices) && 
        length(current_tour_indices) == n_ciudades && 
        all(current_tour_indices > 0) && 
        is.finite(current_cost)) {
      
      # Obtener coordenadas y verificar validez
      ruta_gen_coords <- coords[current_tour_indices, ]
      if(nrow(ruta_gen_coords) == n_ciudades && !any(is.na(ruta_gen_coords$lon))) { 
        
        # Añadir primera ciudad al final para cerrar el ciclo
        ruta_gen_coords <- rbind(ruta_gen_coords, ruta_gen_coords[1, ])
        
        # Crear dataframe de segmentos para esta generación
        segmentos_gen <- data.frame(
          lon_inicio = head(ruta_gen_coords$lon, -1), 
          lat_inicio = head(ruta_gen_coords$lat, -1), 
          lon_fin = tail(ruta_gen_coords$lon, -1), 
          lat_fin = tail(ruta_gen_coords$lat, -1), 
          generacion = gen, 
          costo = current_cost
        )
        
        # Añadir al dataframe general
        all_segments_history_ga <- rbind(all_segments_history_ga, segmentos_gen) 
      } # Fin if coordenadas válidas
    } # Fin if ruta/costo válido
  } # Fin for loop
} # Fin if historial existe

# Verificar si se generaron datos para la animación
if(nrow(all_segments_history_ga) == 0) { 
  print("ADVERTENCIA: No se generaron datos de historial GA para animación.") 
} else { 
  print(paste("Datos de historial GA procesados para", length(unique(all_segments_history_ga$generacion)), "generaciones.")) 
}

## [1] "Datos de historial GA procesados para 200 generaciones."

En la figura 25 se nota como obtiene mejores resultados en especial en el “Best” a medida que el algoritmo avanza en las generaciones.

Creamos la animación para los algortimos geneticos:

# Solo animar si hay datos y mapa base
if(exists("all_segments_history_ga") && nrow(all_segments_history_ga) > 0 && exists("colombia")) {
  
    # 1. Gráfico Base
    mapa_base_historia_ga <- ggplot() + 
      geom_sf(data = colombia, fill = "lightblue", color = "gray50") +
      geom_point(data = coords, aes(x = lon, y = lat), color = "navy", size = 3.0, alpha = 0.6) +
      theme_minimal(base_size = 11) + 
      theme(plot.title = element_text(hjust = 0.5, face="bold"), 
            plot.subtitle = element_text(hjust = 0.5)) +
      coord_sf(xlim = range(coords$lon, na.rm=T) + c(-1.2, 1.2), 
               ylim = range(coords$lat, na.rm=T) + c(-1.2, 1.2), 
               expand = FALSE)

    # 2. Objeto de Animación
    animacion_proceso_ga <- mapa_base_historia_ga +
      geom_segment(
        data = all_segments_history_ga, 
        aes(x = lon_inicio, y = lat_inicio, xend = lon_fin, yend = lat_fin), 
        color = "purple", size = 1.0, alpha = 0.8
      ) + 
      transition_manual(generacion) + 
      labs(
        title = "Proceso GA - Generación {frame}", 
        subtitle = "Mejor ruta encontrada en esta generación", 
        caption = "Figura 26: Proceso de optimización ACO(Fuente: Elaboración propia, 2025)", 
        x = "Longitud", y = "Latitud"
      ) + 
      ease_aes('linear') 
      
    # 3. Parámetros de Renderizado
    fps_proceso_ga <- 7
    total_frames_proceso_ga <- length(unique(all_segments_history_ga$generacion))
    
    print(paste("Renderizando animación GA:", total_frames_proceso_ga, "frames a", fps_proceso_ga, "fps."))
    
    # 4. Renderizar (con tryCatch)
    anim_render_proceso_ga <- tryCatch({
      animate(
        animacion_proceso_ga, 
        nframes = total_frames_proceso_ga, 
        fps = fps_proceso_ga, 
        width = 800, 
        height = 700, 
        renderer = gifski_renderer()
      )
    }, error=function(e){
        print(paste("Error renderizando animación GA:", e$message))
        return(NULL)
    })
    
    # 5. Visualizar y Guardar (si renderización fue exitosa)
    if (!is.null(anim_render_proceso_ga)) { 
      # Mostrar en dispositivo gráfico / Rmd
      print(anim_render_proceso_ga)
      
      # Guardar
      tryCatch({
          anim_save("proceso_optimizacion_ga.gif", animation = anim_render_proceso_ga)
          print("Animación Proceso GA guardada.") 
      }, error=function(e) {
          print(paste("Error guardando animación GA:", e$message))
      })
    } else { 
      print("Fallo al renderizar animación GA.") 
    }
    
} else { 
  print("No hay datos/mapa para generar animación GA.") 
}

## [1] "Renderizando animación GA: 200 frames a 7 fps."

## Warning in device(files[i], width = 800, height = 700, units = "in", res =
## 192): cairo error 'invalid value (typically too big) for the size of the input
## (surface, pattern, etc.)'

## [1] "Error renderizando animación GA: unable to start device 'png'"
## [1] "Fallo al renderizar animación GA."

La animación resultante (Figura 26) ilustra cómo la mejor ruta encontrada por el algoritmo AG evoluciona a lo largo de las iteraciones, tendiendo a estabilizarse a medida que el algoritmo converge hacia una solución de bajo costo.

knitr::include_graphics("proceso_optimizacion_ga.gif")

Comparando ambos algoritmos:

costo_aco <-
  if (!is.null(resultado_aco) && is.finite(resultado_aco$best_tour_cost)) {
    resultado_aco$best_tour_cost
  } else {
    Inf # Asignar Infinito si no hay resultado válido
  }
indices_aco <-
  if (!is.null(resultado_aco) && is.finite(resultado_aco$best_tour_cost)) {
    resultado_aco$best_tour_indices
  } else {
    NULL # Asignar NULL si no hay resultado válido
  }

costo_ga <- if (exists("best_tour_cost_ga") && is.finite(best_tour_cost_ga)) best_tour_cost_ga else Inf
indices_ga <- if (exists("best_tour_indices_ga") && is.finite(best_tour_cost_ga)) best_tour_indices_ga else NULL

if (is.infinite(costo_aco) && is.infinite(costo_ga)) {
  print("Error: Ninguno de los algoritmos encontró una solución válida.")
  final_best_cost <- Inf
  final_best_tour_indices <- NULL
  algoritmo_ganador <- "Ninguno"
} else {
  if (costo_ga < costo_aco) {
    algoritmo_ganador <- "Algoritmo Genético (GA)"
    final_best_cost <- costo_ga
    final_best_tour_indices <- indices_ga
    print(paste("GA encontró una mejor solución (Costo:", format(final_best_cost, scientific = FALSE, big.mark=","), ")"))
    print(paste("  (Costo ACO:", format(costo_aco, scientific = FALSE, big.mark=","), ")"))
  } else if (costo_aco < costo_ga) {
    algoritmo_ganador <- "Colonia de Hormigas (ACO)"
    final_best_cost <- costo_aco
    final_best_tour_indices <- indices_aco
    print(paste("ACO encontró una mejor solución (Costo:", format(final_best_cost, scientific = FALSE, big.mark=","), ")"))
    print(paste("  (Costo GA:", format(costo_ga, scientific = FALSE, big.mark=","), ")"))
  } else { # Costos iguales
    algoritmo_ganador <- "Ambos (costo igual)"
    final_best_cost <- costo_ga # O costo_aco, son iguales
    # Podemos elegir arbitrariamente GA en caso de empate, o mantener ambas si quisiéramos
    final_best_tour_indices <- indices_ga 
    print(paste("Ambos algoritmos encontraron soluciones con el mismo costo:", format(final_best_cost, scientific = FALSE, big.mark=",")))
  }
  
  # Imprimir la mejor ruta encontrada
  if (!is.null(final_best_tour_indices)) {
      final_best_tour_cities <- ciudades_nombres[final_best_tour_indices]
      print(paste("Mejor ruta encontrada por", algoritmo_ganador, ":"))
      print(paste(final_best_tour_cities, collapse = " -> "))
  }
}

## [1] "GA encontró una mejor solución (Costo: 2,403,978 )"
## [1] "  (Costo ACO: 2,860,689 )"
## [1] "Mejor ruta encontrada por Algoritmo Genético (GA) :"
## [1] "Medellin -> Manizales -> Pereira -> Cali -> Pasto -> Neiva -> Ibague -> Bogota -> Bucaramanga -> Cucuta -> Santa Marta -> Barranquilla -> Cartagena"

En conclusión, el GA superó al ACO en ≈ 16 % de ahorro de costo, lo que indica que, con la configuración actual, el GA exploró el espacio de soluciones de forma más efectiva para este problema de rutas entre 13 ciudades.

La mejor ruta hallada por el GA cubre primero el suroccidente (Pasto → Cali → Pereira → Manizales → Medellín) y luego “salta” a la costa Caribe (Cartagena → Barranquilla → Santa Marta), para finalmente recorrer el nororiente y centro (Cúcuta → Bucaramanga → Bogotá → Ibagué → Neiva).
Ese orden evita retornos largos y, al juntar tramos geográficamente cercanos, reduce tanto distancias como tiempo de conducción.

En contraste, el ACO depende fuertemente de los parámetros de feromonas (α, β, tasa de evaporación ). Si la evaporación es alta o la influencia heurística (β) es baja, las hormigas pueden dispersarse demasiado y tardar más en reforzar la mejor ruta. Con los valores actuales, el ACO siguió mejorando pero no alcanzó la calidad que el GA obtuvo antes de que se cumpliera su criterio de parada.

Curso: Redes Neuronales Artificiales y Algoritmos Bioinspirados

Semestre: 2025 - 1

Trabajo #1: Optimización Heurística

Equipo #2:

Tabla de contenidos

Parte 1: Optimización Numérica

Funciones de Prueba

Optimización

Descenso por gradiente

Rastrigin

Algoritmos evolutivos

Rastrigin

Rosenbrock

Enjambres de partículas

Rastrigin

Rosenbrock

Evolución diferencial

Rastrigin

Rosenborck

Generación de los GIF

Descenso por gradiente

Algoritmos evolutivos

Enjambres de partículas

Evolución diferencial

Heuristícos VS Descenso por gradiente

Rastrigin

Rosenbrock

Discusión

Parte 2: Optimización Combinatoria

Descripción del problema del vendedor

Algoritmos utilizados

Código en R para la optimización

Algoritmo de Optimización por Colonia de Hormigas (ACO)

Bibliografía

Reporte de Contribución Individual