Trabajo 1: Optimización Heurística

REDES NEURONALES Y ALGORITMOS BIOINSPIRADOS

Presentado por:

Juan Jose Correa Hurtado <jjcorreahu@unal.edu.co>,

Jacobo Ochoa Ramirez <jochoara@unal.edu.co>,

Hailer Serna Hernandez <hsernah@unal.edu.co>

Profesor: Juan David Ospina Arango

Monitor: Andrés Mauricio Zapata Rincón

Universidad Nacional de Colombia Facultad de Minas Ingeniería de Sistemas e Informática

02 de mayo de 2025

Importación de Librerias

library(GA)
library(ggplot2)
library(reshape2)
library(plot3D)
library(animation)
library(DEoptim)
library(pso)
library(dplyr)
library(gifski)
library(gganimate)
library(knitr)
library(magick)

OPTIMIZACIÓN NUMÉRICA

Para evaluar el desempeño de distintos métodos de optimización numérica, se seleccionaron dos funciones de prueba ampliamente utilizadas en la literatura: la función de Rosenbrock y la función de Schwefel. Estas funciones presentan características contrastantes que permiten observar cómo los algoritmos se comportan en contextos diferentes: con mínimos globales difíciles de alcanzar o con múltiples mínimos locales que pueden confundir la convergencia.

Función de Rosenbrock

La función de Rosenbrock es una función unimodal cuyo mínimo global se encuentra en un valle angosto y curvado de forma parabólica. Aunque este valle es fácil de localizar visualmente, alcanzar el mínimo representa un reto para muchos algoritmos debido a la geometría del terreno de optimización.

“La función es unimodal, y el mínimo global se encuentra en un valle parabólico angosto. Sin embargo, aunque este valle es fácil de encontrar, la convergencia al mínimo es difícil” (Surjanovic & Bingham, s.f.) [traducción propia].

Esto la convierte en un excelente caso para probar algoritmos que deben navegar por terrenos suaves pero estrechos, donde los gradientes pueden apuntar en direcciones poco útiles si no se ajustan bien los parámetros de convergencia.

Para \(d\) dimensiones la función de Rosenbrock está dada por la forma:

\[ f(\mathbf{x}) = \sum_{i=1}^{d-1} \left[100(x_{i+1} - x_i^2)^2 + (x_i-1)^2\right] \]

El siguiente código implementa la función:

funcion_rosenbrock <- function(vector_x){
  
  # Evaluación de la función de Rosenbrock para un vector de entrada vector_x
  
  # Se expresa de esta manera para prevenir el uso de ciclos y 
  # sacar mayor provecho al paradigma de programación funcional
  
  d <- length(vector_x)   # Obtener la dimensión del vector
  xi <- vector_x[1:(d-1)] # x_i
  xnext <- vector_x[2:d]  # x_(i+1)
  
  # f(x) = sum_{i=1}^{d-1} [100 * (x_{i+1} - x_i^2)^2 + (1 - x_i)^2]
  result <- (sum(100*(xnext-xi^2)^2+(xi-1)^2))
  return(result)
}

Representación 3D de la función:

Función de Schwefel

En contraste, la función de Schwefel presenta una superficie altamente compleja, con múltiples mínimos locales dispersos en todo el dominio. Esta característica desafía a los métodos de optimización a evitar caer en soluciones subóptimas y a explorar eficientemente el espacio de búsqueda.

“La función de Schwefel es compleja, con muchos mínimos locales” (Surjanovic & Bingham, s.f.) [traducción propia].

Este tipo de función es ideal para evaluar algoritmos en escenarios donde la exploración global es tan importante como la exploración local.

Para \(d\) dimensiones la función de Schwefel está dada por la forma:

\[ f(\mathbf{x}) = 418.9829 \cdot d - \sum_{i=1}^{d} x_i \sin\left(\sqrt{|x_i|}\right) \]

El siguiente código implementa la función:

funcion_schwefel <- function(vector_x){
  
  # Evaluación de la función Schwefel dado un vector de entrada vector_x 
  
  d <- length(vector_x) # Obtener la dimensión del vector
  
  
  # f(x) = 418.9829 * d - sum(x_i * sin(sqrt(|x_i|)))
  # Este valor tiene su mínimo global en x_i ≈ 420.9687 para todo i
  
  sum <- sum(vector_x*sin(sqrt(abs(vector_x))))
  result <- 418.9829*d - sum
  
  return(result)
}

Representación 3D de la función:

Optimización por Descenso del Gradiente

Para poder optimizar ambas funciones en dos y tres dimensiones se necesita un optimizador que haga uso del gradiente numérico.

Inicialmente se debe obtener la derivada parcial de la función respecto a alguna de sus componentes.

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

Se evalua el gradiente de la función al obtener cada una de las derivadas parciales de \(f\) en \(x\)

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

EL siguiente código hace uso del gradiente númerico para calcular el descenso del gradiente númerico con una tasa de aprendizaje variable:

optimizador_multivariado <- function(fun, x0, eta0=0.1, decay=0.01, max_eval=100, h=0.01, tol = 1e-5) {


  x <- matrix(NA, nrow=max_eval, ncol=length(x0))  # Matriz para registrar los puntos iterativos
  f_values <- numeric(max_eval)       # Vector para registrar los valores de la función objetivo
  
  
  x[1, ] <- x0                # asignar a la primer columna de la matriz x el vector x0
  f_values[1] <- fun(x0)      # valor inicial de la función
  
  
  for (i in 2:max_eval) {
    
    # Calcula el gradiente numérico en la iteración anterior
    grad <- num_grad(x[i-1, ], fun, h)   
    
    # Tasa de aprendizaje variable (disminuye con el tiempo)
    eta <- eta0 / (1 + decay * (i - 1)) 
    
    # Paso de descenso: actualiza x restando el gradiente
    x[i, ] <- x[i-1, ] - eta * grad      
    
    # Calcular valor de la función en el nuevo punto
    f_values[i] <- fun(x[i, ])
    
    # Magnitud del cambio entre iteraciones
    cambio <- sqrt(sum((x[i, ] - x[i-1, ])^2))     
    
    # Si el cambio es muy pequeño, se asume convergencia
    if (cambio < tol) {                
      break
    }
  }
  
  
  result <- list(
    puntos = x[1:i, ],          # Puntos recorridos
    valores = f_values[1:i],    # Valores de la función objetivo
    iteraciones = i,            # Número total de iteraciones
    solucion = x[i, ],          # Solución final encontrada
    valor_optimo = f_values[i]  # Valor óptimo encontrado
  )
  
  return(result)  
}

Condición Inicial Aleatoria

EL siguiente código genera un vector de \(n\) dimensiones para ser evaluado en el optimizador multivariado.

vector_aleatorio <- function(n,limInf, limSup, seed = NULL){
  if (!is.null(seed)){
    set.seed(seed) # semilla fija solo si fue proporcionada 
  }
  
  return(runif(n, min = limInf, max = limSup)) # Vector aleatorio en el rango proporcionado
}

Optimización de Función Rosenbrock

Dos Dimensiones:

x0_2d_rosenbrock <- vector_aleatorio(2, -2.048 , 2.048, seed = 28)

res_rosenbrock_2d <- optimizador_multivariado(funcion_rosenbrock, 
                                              x0 = x0_2d_rosenbrock, 
                                              eta0 = 0.001,
                                              decay = 0.2)

cat("Solución: ",res_rosenbrock_2d$solucion, "\nValor optimo: ", res_rosenbrock_2d$valor_optimo)
## Solución:  0.08757869 -0.01458303 
## Valor optimo:  0.8820325
sol <- res_rosenbrock_2d$puntos

n_length <- 100
x1 <- seq(-2, 2, length.out = n_length)
x2 <- seq(-2, 2, length.out = n_length)
X <- expand.grid(x1, x2)

# Evaluar la función en toda la grilla
z <- apply(X, 1, function(row) funcion_rosenbrock(c(row[1], row[2])))
Z <- matrix(z, ncol = n_length, nrow = n_length)

# Crear gráfico de curvas de nivel
contour(x1, x2, Z,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Ejemplo: Curvas de nivel - Función de Rosenbrock",
        sub = "Camino del optimizador multivariado",
        las = 1, drawlabels = TRUE)

#leyendas
legend(-1.5, 2, legend = c("Punto de inicio", "Pasos"), fill = c("blue", "red"))
#Pasos
lines(sol[,1], sol[,2], type = "b", pch = 19, col = "red", lwd = 2)
#Punto de incio
points(sol[1,1], sol[1,2], col = "blue", pch = 19, cex = 1.8)

Tres Dimensiones:

x0_3d_rosenbrock <- vector_aleatorio(3, -2.048 , 2.048, seed = 45)

res_rosenbrock_3d <- optimizador_multivariado(funcion_rosenbrock, 
                                              x0 = x0_3d_rosenbrock, 
                                              eta0 = 0.001,
                                              decay = 0.2)

cat("Solución: ",res_rosenbrock_3d$solucion, "\nValor optimo: ", res_rosenbrock_3d$valor_optimo)
## Solución:  0.2708469 0.06569966 -0.04255669 
## Valor optimo:  1.630156

Optimización de Función de Schwefel

Dos Dimensiones:

x0_2d_schwefel <- vector_aleatorio(2, -500 , 500, seed = 69)

res_schwefel_2d <- optimizador_multivariado(funcion_schwefel, 
                                              x0 = x0_2d_schwefel, 
                                              eta0 = 0.3,
                                              decay = 0.01)

cat("Vector origen: ", x0_2d_schwefel,"\nSolucion: ",res_schwefel_2d$solucion, "\nValor optimo: ", res_schwefel_2d$valor_optimo)
## Vector origen:  30.75401 268.8077 
## Solucion:  65.07363 204.3101 
## Valor optimo:  572.5487
sol <- res_schwefel_2d$puntos

n_length <- 100
x1 <- seq(-500, 500, length.out = n_length)
x2 <- seq(-500, 500, length.out = n_length)
X <- expand.grid(x1, x2)

# Evaluar la función en toda la grilla
z <- apply(X, 1, function(row) funcion_schwefel(c(row[1], row[2])))
Z <- matrix(z, ncol = n_length, nrow = n_length)

# Crear gráfico de curvas de nivel
contour(x1, x2, Z,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Curvas de nivel - Funcion de schwefel",
        sub = "Camino del optimizador multivariado",
        las = 1, drawlabels = TRUE)

#leyendas
legend(-400, -200, legend = c("Punto de inicio", "Pasos"), fill = c("blue", "red"))
#Pasos
lines(sol[,1], sol[,2], type = "b", pch = 19, col = "red", lwd = 2)
#Punto de incio
points(sol[1,1], sol[1,2], col = "blue", pch = 19, cex = 1.8)

sol <- res_schwefel_2d$puntos

n_length <- 100
x1 <- seq(0, 100, length.out = n_length)
x2 <- seq(180, 300, length.out = n_length)
X <- expand.grid(x1, x2)

# Evaluar la función en toda la grilla
z <- apply(X, 1, function(row) funcion_schwefel(c(row[1], row[2])))
Z <- matrix(z, ncol = n_length, nrow = n_length)

# Crear gráfico de curvas de nivel
contour(x1, x2, Z,
        xlab = expression(x[1]), ylab = expression(x[2]),
        main = "Curvas de nivel con zoom - Funcion de schwefel",
        sub = "Camino del optimizador multivariado",
        las = 1, drawlabels = TRUE)

#leyendas
legend(80, 290, legend = c("Punto de inicio", "Pasos"), fill = c("blue", "red"))
#Pasos
lines(sol[,1], sol[,2], type = "b", pch = 19, col = "red", lwd = 2)
#Punto de incio
points(sol[1,1], sol[1,2], col = "blue", pch = 19, cex = 1.8)

Tres Dimensiones:

x0_3d_schwefel <- vector_aleatorio(3, -2.048 , 2.048, seed = 45)

res_schwefel_3d <- optimizador_multivariado(funcion_schwefel, 
                                              x0 = x0_3d_schwefel, 
                                              eta0 = 0.3,
                                              decay = 0.01)

cat("Solución: ",res_schwefel_3d$solucion, "\nValor optimo: ", res_schwefel_3d$valor_optimo)
## Solución:  5.238147 5.236521 5.236186 
## Valor optimo:  1245.113

Optimización Heurística

Ante las limitaciones del algoritmo de optimización por descenso del gradiente se plantea el uso de otros tres algoritmos que pueden tener un mejor desempeño.Estos son:

  • Algoritmos Evolutivos.
  • Optimización de Particulas.
  • Evolución Diferencial.

Algoritmos Evolutivos

Los algoritmos evolutivos son técnicas de optimización inspiradas en los principios de la selección natural y genética. Operan sobre una población de soluciones que evoluciona a través de generaciones mediante operadores como selección, cruzamiento y mutación. Cada individuo representa una solución candidata al problema, y su adecuación se evalúa mediante una función objetivo. A lo largo de las iteraciones, las soluciones de mayor aptitud tienden a prevalecer, promoviendo la exploración y explotación del espacio de búsqueda.

Mediante el uso de la librería 'GA' se va implementar un algoritmo genético para optimizar ambas funciones objetivo.

Para evaluar el comportamiento del algoritmo post-evaluación se va a definir una función que registra la población de cada generación junto con el mejor individuo de la población.

make_postfit <- function(pop_store_name, best_store_name = NULL) {
  # Clear previous population variable
  assign(pop_store_name, list(), envir = globalenv())

  # Clear best individual variable if provided
  if (!is.null(best_store_name)) {
    assign(best_store_name, list(), envir = globalenv())
  }

  function(object, ...) {
    pop <- object@population
    fitness <- object@fitness

    # Append current population
    pop_list <- get(pop_store_name, envir = globalenv())
    pop_list <- append(pop_list, list(pop))
    assign(pop_store_name, pop_list, envir = globalenv())

    # If best_store_name is provided, track best individual
    if (!is.null(best_store_name)) {
      best_index <- which.max(fitness)  # For maximization
      best_individual <- pop[best_index, , drop = FALSE]

      best_list <- get(best_store_name, envir = globalenv())
      best_list <- append(best_list, list(best_individual))
      assign(best_store_name, best_list, envir = globalenv())
    }

    return(object)
  }
}

Optimización Mediante Algoritmo Genético de Función de Rosenbrock

Dos Dimensiones:

Se procede a realizar la optimización en dos dimensiones de la función Rosenbrock.Se va a hacer uso de algoritmos bioinspirados como selección, mutación y cruzamiento, cada uno con sus respectivos parametros, adecuados para obtener mejores resultados.

postfit_rosenbrock_2d <- make_postfit("pop_rosenbrock", "best_rosenbrock")

TYPE = "real-valued"
popSize = 50

GA_rosenbrock <- ga(type = TYPE,
                    fitness = function(x){-funcion_rosenbrock(x)},
                    lower = c(-2.048 , -2.048), upper = c(2.048 , 2.048),
                    maxiter = 1000,
                    run = 100,
                    population = gaControl(TYPE)$population,
                    selection = gaControl(TYPE)$selection,
                    crossover = gaControl(TYPE)$crossover, 
                    mutation = gaControl(TYPE)$mutation,
                    popSize = popSize,
                    pcrossover = 0.88, 
                    pmutation = 0.1, 
                    elitism = base::max(1, round(popSize*0.05)), 
                    seed = 555,
                    postFitness = postfit_rosenbrock_2d
                  
                    )
summary(GA_rosenbrock)
## ── Genetic Algorithm ─────────────────── 
## 
## GA settings: 
## Type                  =  real-valued 
## Population size       =  50 
## Number of generations =  1000 
## Elitism               =  2 
## Crossover probability =  0.88 
## Mutation probability  =  0.1 
## Search domain = 
##           x1     x2
## lower -2.048 -2.048
## upper  2.048  2.048
## 
## GA results: 
## Iterations             = 256 
## Fitness function value = -0.003070611 
## Solution = 
##             x1        x2
## [1,] 0.9449736 0.8923217
plot(GA_rosenbrock)

x1 <- x2 <- seq(-5.12, 5.12, by = 0.1)

adapt_rosenbrock <- function(x1, x2) {
  mapply(function(a, b) funcion_rosenbrock(c(a, b)), x1, x2)
}
f <- outer(x1, x2, adapt_rosenbrock)
iter_to_show = c(1,60,120, 180, 240, 256)
par(mfrow = c(3,2), mar = c(2,2,2,1),
    mgp = c(1, 0.4, 0), tck = -.01)
for(i in seq(iter_to_show))
{
  contour(x1, x2, f, drawlabels = FALSE, col = "grey50")
  title(paste("Iter =", iter_to_show[i]))
  points(pop_rosenbrock[[iter_to_show[i]]], pch = 20, col = "forestgreen")
}

Tres Dimensiones:

Ahora se va a realizar la optimización para tres dimensiones de la función de Rosenbrock.

postfit_rosenbrock_3d <- make_postfit("pop_rosenbrock_3d", "best_rosenbrock_3d")

TYPE = "real-valued"
popSize = 200

GA_rosenbrock_3d <- ga(type = TYPE,
                    fitness = function(x){-funcion_rosenbrock(x)},
                    lower = c(-2.048 , -2.048, -2.048), upper = c(2.048 , 2.048, 2.048),
                    maxiter = 2000,
                    run = 200,
                    population = gaControl(TYPE)$population,
                    selection = gaControl(TYPE)$selection,
                    crossover = gaControl(TYPE)$crossover, 
                    mutation = gaControl(TYPE)$mutation,
                    popSize = popSize,
                    pcrossover = 0.64, 
                    pmutation = 0.05, 
                    elitism = base::max(1, round(popSize*0.13)), 
                    seed = 555,
                    postFitness = postfit_rosenbrock_3d
                  
                    )
summary(GA_rosenbrock_3d)
## ── Genetic Algorithm ─────────────────── 
## 
## GA settings: 
## Type                  =  real-valued 
## Population size       =  200 
## Number of generations =  2000 
## Elitism               =  26 
## Crossover probability =  0.64 
## Mutation probability  =  0.05 
## Search domain = 
##           x1     x2     x3
## lower -2.048 -2.048 -2.048
## upper  2.048  2.048  2.048
## 
## GA results: 
## Iterations             = 2000 
## Fitness function value = -0.1377346 
## Solution = 
##            x1        x2       x3
## [1,] 0.821034 0.6753304 0.454874

Optimización Mediante Algoritmo Genético de Función de Schwefel

Dos Dimensiones:

Se procede a realizar la optimización en dos dimensiones de la función de Schwefel.

postfit_schwefel_2d <- make_postfit("pop_schwefel_2d", "best_schwefel_2d")

TYPE = "real-valued"
popSize = 50

GA_schwefel <- ga(type = TYPE,
                    fitness = function(x){-funcion_schwefel(x)},
                    lower = c(-500 , -500), upper = c(500 , 500),
                    maxiter = 1000,
                    run = 100,
                    population = gaControl(TYPE)$population,
                    selection = gaControl(TYPE)$selection,
                    crossover = gaControl(TYPE)$crossover, 
                    mutation = gaControl(TYPE)$mutation,
                    popSize = popSize,
                    pcrossover = 0.88, 
                    pmutation = 0.1, 
                    elitism = base::max(1, round(popSize*0.05)), 
                    seed = 555,
                    postFitness = postfit_schwefel_2d
                  
                    )
summary(GA_schwefel)
## ── Genetic Algorithm ─────────────────── 
## 
## GA settings: 
## Type                  =  real-valued 
## Population size       =  50 
## Number of generations =  1000 
## Elitism               =  2 
## Crossover probability =  0.88 
## Mutation probability  =  0.1 
## Search domain = 
##         x1   x2
## lower -500 -500
## upper  500  500
## 
## GA results: 
## Iterations             = 241 
## Fitness function value = -0.000586059 
## Solution = 
##            x1       x2
## [1,] 420.9069 420.9438

Tres Dimensiones

postfit_schwefel_3d <- make_postfit("pop_schwefel_3d", "best_schwefel_3d")

TYPE = "real-valued"
popSize = 200

GA_schwefel_3d <- ga(type = TYPE,
                  fitness = function(x) {-funcion_schwefel(x)},
                  lower = c(-500, -500, -500),
                  upper = c(500, 500, 500),
                  maxiter = 1000,
                  run = 100,
                  population = gaControl(TYPE)$population,
                  selection = gaControl(TYPE)$selection,
                  crossover = gaControl(TYPE)$crossover, 
                  mutation = gaControl(TYPE)$mutation,
                  popSize = popSize,
                  pcrossover = 0.64, 
                  pmutation = 0.05, 
                  elitism = base::max(1, round(popSize * 0.13)), 
                  seed = 555,
                  postFitness = postfit_schwefel_3d
)
summary(GA_schwefel_3d)
## ── Genetic Algorithm ─────────────────── 
## 
## GA settings: 
## Type                  =  real-valued 
## Population size       =  200 
## Number of generations =  1000 
## Elitism               =  26 
## Crossover probability =  0.64 
## Mutation probability  =  0.05 
## Search domain = 
##         x1   x2   x3
## lower -500 -500 -500
## upper  500  500  500
## 
## GA results: 
## Iterations             = 334 
## Fitness function value = -3.81827e-05 
## Solutions = 
##             x1       x2       x3
## [1,]  420.9687 420.9687 420.9687
## [2,]  420.9687 420.9687 420.9687
## [3,]  420.9687 420.9687 420.9687
## [4,]  420.9687 420.9687 420.9687
## [5,]  420.9687 420.9687 420.9687
## [6,]  420.9687 420.9687 420.9687
## [7,]  420.9687 420.9687 420.9687
## [8,]  420.9687 420.9687 420.9687
## [9,]  420.9687 420.9687 420.9687
## [10,] 420.9687 420.9687 420.9687
##  ...                            
## [26,] 420.9687 420.9687 420.9687
## [27,] 420.9687 420.9687 420.9687

Optimización de Particulas

La optimización por enjambre de partículas es un algoritmo bioinspirado basado en el comportamiento colectivo observado en sistemas naturales como bandadas de aves o bancos de peces. Cada partícula representa una solución potencial que “vuela” en el espacio de búsqueda ajustando su posición según su experiencia personal y la de sus vecinas. Las actualizaciones se rigen por la mejor posición individual y global conocidas, lo que permite un equilibrio entre la búsqueda local y global de óptimos.

Mediante el uso de la librería 'pso' se va implementar el algoritmo para optimizar las dos funciones en cuestion.

# Particle Swarm Optimization Function
pso_optimization <- function(
  func, # La función a optimizar
  dim,   # Dimensión
  n_particles = 30, # Número de partículas en el enjambre
  max_iter = 100,  # Número máximo de iteraciones
  w = 0.7,         # Peso de la inercia
  c1 = 1.4,        # Coeficiente de cognitivo 
  c2 = 1.4,        # Coeficiente social
  lower_bound,     # Límite inferior del espacio de búsqueda
  upper_bound,     # Límite superior del espacio de búsqueda
  target_fitness = -Inf # Valor de aptitud física objetivo para detenerse temprano.
) {

  # Inicializar partículas y velocidades 
  particles <- matrix(runif(n_particles * dim, min = lower_bound, max = upper_bound),
                      nrow = n_particles, ncol = dim)
  velocities <- matrix(runif(n_particles * dim, min = -0.1 * (upper_bound - lower_bound),
                             max = 0.1 * (upper_bound - lower_bound)),
                       nrow = n_particles, ncol = dim)
  personal_best_positions <- particles
  personal_best_fitnesses <- apply(particles, 1, func) # Use apply for row-wise fitness evaluation
  global_best_position <- particles[which.min(personal_best_fitnesses), ]
  global_best_fitness <- min(personal_best_fitnesses)

  # Almacenar posiciones de partículas para cada iteración para graficar
  particle_positions_history <- list()
  fitness_history <- c() # track the best fitness over iterations

  # Bucle de optimización
  for (iter in 1:max_iter) {
    # Update velocities and positions
    r1 <- matrix(runif(n_particles * dim), nrow = n_particles, ncol = dim)
    r2 <- matrix(runif(n_particles * dim), nrow = n_particles, ncol = dim)

    velocities <- w * velocities +
      c1 * r1 * (personal_best_positions - particles) +
      c2 * r2 * (matrix(rep(global_best_position, n_particles),
                       nrow = n_particles, byrow = TRUE) - particles)
    particles <- particles + velocities

    # Recortar partículas a los límites del espacio de búsqueda
    particles <- pmax(particles, lower_bound)
    particles <- pmin(particles, upper_bound)

    # Evaluar la idoneidad de los nuevos puestos
    current_fitnesses <- apply(particles, 1, func) # Use apply here too

    # Actualizar mejores posiciones y estados físicos personales
    improved_mask <- current_fitnesses < personal_best_fitnesses
    personal_best_positions[improved_mask, ] <- particles[improved_mask, ]
    personal_best_fitnesses[improved_mask] <- current_fitnesses[improved_mask]

    # Actualizar la mejor posición y condición física a nivel mundial
    if (min(current_fitnesses) < global_best_fitness) {
      global_best_position <- particles[which.min(current_fitnesses), ]
      global_best_fitness <- min(current_fitnesses)
    }

    # Almacenar posiciones de partículas para esta iteración
    particle_positions_history[[iter]] <- particles
    fitness_history <- c(fitness_history, global_best_fitness)

    # Compruebe si hay parada anticipada
    if (global_best_fitness <= target_fitness) {
      break
    }
  }

  return(list(
    particles = particles,
    global_best_position = global_best_position,
    global_best_fitness = global_best_fitness,
    particle_positions_history = particle_positions_history, # Return history
    fitness_history = fitness_history,
    iterations = iter
  ))
}
# Parámetros
dim <- 2
n_particles <- 50
max_iter <- 10
lower_bound <- -5
upper_bound <- 5

# Ejecute PSO en la función Rosenbrock
result_rosenbrock <- pso_optimization(funcion_rosenbrock, dim, n_particles, max_iter,
                                     lower_bound = lower_bound, upper_bound = upper_bound)

# Ejecute PSO en la función Schwefel
result_schwefel <- pso_optimization(funcion_schwefel, dim, n_particles, max_iter,
                                    lower_bound = -500, upper_bound = 500, 
                                    target_fitness = 0) # Example target fitness

Evolución Diferencial

La evolución diferencial es una técnica estocástica de optimización continua que opera sobre una población de vectores mediante mecanismos de mutación diferencial, recombinación y selección. Su característica distintiva es la generación de nuevos candidatos mediante combinaciones lineales de soluciones existentes, lo que facilita una exploración eficiente del espacio de búsqueda. DE ha demostrado ser particularmente eficaz en problemas no lineales, no diferenciables y multimodales.

# OPTIMIZACIÓN de la función
differential_evolution <- function(fn, bounds, max_iter = 100, pop_size = 30, F = 0.8, CR = 0.9) {
  D <- length(bounds[[1]])
  pop <- matrix(runif(pop_size * D, bounds[[1]], bounds[[2]]), nrow = pop_size, byrow = TRUE)
  fitness <- apply(pop, 1, fn)
  history <- list(pop)  # guardar posiciones

  for (iter in 1:max_iter) {
    for (i in 1:pop_size) {
      idxs <- sample(setdiff(1:pop_size, i), 3)
      x1 <- pop[idxs[1],]
      x2 <- pop[idxs[2],]
      x3 <- pop[idxs[3],]
      mutant <- x1 + F * (x2 - x3)
      mutant <- pmin(pmax(mutant, bounds[[1]]), bounds[[2]])  # clip

      crossover <- runif(D) < CR
      if (!any(crossover)) crossover[sample(1:D, 1)] <- TRUE
      trial <- ifelse(crossover, mutant, pop[i,])

      f_trial <- fn(trial)
      if (f_trial < fitness[i]) {
        pop[i,] <- trial
        fitness[i] <- f_trial
      }
    }
    history[[iter + 1]] <- pop
  }
  list(best = pop[which.min(fitness),], value = min(fitness), history = history)
}

# Visualización de la optimización para 2D
plot_history_gif <- function(history, fn, bounds, filename = "de_optimization.gif") {
  frames <- list()
  grid_points <- expand.grid(
    x = seq(bounds[[1]][1], bounds[[2]][1], length.out = 100),
    y = seq(bounds[[1]][2], bounds[[2]][2], length.out = 100)
  )
  grid_points$z <- apply(grid_points, 1, fn)

  for (i in seq_along(history)) {
    df <- as.data.frame(history[[i]])
    colnames(df) <- c("x", "y")
    p <- ggplot() +
      geom_raster(data = grid_points, aes(x, y, fill = z), alpha = 0.5) +
      scale_fill_viridis_c() +
      geom_point(data = df, aes(x, y), color = "red", size = 2) +
      ggtitle(sprintf("Iteración %d", i)) +
      xlim(bounds[[1]][1], bounds[[2]][1]) +
      ylim(bounds[[1]][2], bounds[[2]][2]) +
      theme_minimal()
    img <- image_graph(600, 600, res = 96)
    print(p)
    frames[[i]] <- image_capture()
    dev.off()
  }

  animation <- image_animate(image_join(frames), fps = 5)
  image_write(animation,filename)
}
#USO DE LAS FUNCIONES Y GENERACIÓN DE GIFS
set.seed(1983)

# Elegir función y dominio

bounds_rs <- list(c(-2, -2), c(2, 2))  # para Rosenbrock
bounds_sc <- list(c(-500, -500), c(500, 500))  # para Schwefel

result_rs <- differential_evolution(funcion_rosenbrock, bounds_rs, max_iter = 50, pop_size = 30)


result_sc <- differential_evolution(funcion_schwefel, bounds_sc, max_iter = 50, pop_size = 30)
plot_history_gif(result_sc$history, funcion_schwefel, bounds_sc, filename = "optimizacion_DE_sc.gif")
plot_history_gif(result_rs$history, funcion_rosenbrock, bounds_rs, filename = "optimizacion_DE_rs.gif")
cat("Mejor solución encontrada Rosenbrock:\n")
## Mejor solución encontrada Rosenbrock:
print(result_rs$best)
## [1] 0.9999241 0.9999315
cat("Mejor solución encontrada Schwefel:\n")
## Mejor solución encontrada Schwefel:
print(result_sc$best)
## [1] 421.1128 421.1281
cat("Valor óptimo Rosenbrock:\n")
## Valor óptimo Rosenbrock:
print(result_rs$value)
## [1] 7.002326e-07
cat("Valor óptimo Schwefel:\n")
## Valor óptimo Schwefel:
print(result_sc$value)
## [1] 0.005848507

Comparación de GIFs

Rosenbrock PSO

##### Schwefel PSO ##### Schwefel GA 3D ##### Schwefel GA 2D

knitr::include_graphics("GA_Schwefel.gif")

Rosenbrock GA 3D

Análisis de Animaciones

Para funciones suaves y bien comportadas como Rosenbrock, el descenso por gradiente es eficiente y preciso, siempre que se parta de una buena condición inicial. Mientras que para funciones altamente multimodales o con formas complejas como Schwefel, los métodos heurísticos son superiores en la exploración global y probabilidad de encontrar el mínimo global.

En la práctica, una estrategia mixta (heurístico global + refinamiento con gradiente) puede ofrecer lo mejor de ambos mundos. Dependiendo del contexto del problema, es probable que un método sea mejor que otro, los métodos heurísticos tienen más variabilidad y puede que tome más tiempo encontrar una solución que cumpla con los requisitos establecidos pero en caso de tenerse clara las condiciones iniciales que tiene el problema quizás sea mejor un descenso por gradiente, esto también puede ser algo negativo, el descenso por gradiente es sensible a estas condiciones iniciales. En conclusión, al momento de elegir un método es importante tener un entendimiento de tanto las opciones disponibles a aplicar como la complejidad del problema, es importante usar la herramienta adecuada, ni la más simple ni la más compleja.

# OPTIMIZACIÓN COMBINATORIA

En este documento se resuelve el problema del vendedor viajero visitando las 13 ciudades principales de Colombia. Para encontrar la mejor ruta se utilizan dos métodos de optimización: un Algoritmo Genético y un Algoritmo de Colonias de Hormigas. Se construye una matriz de costos en pesos colombianos, donde todo se traduce a valores monetarios: el costo del combustible según la distancia recorrida, el valor de los peajes y el tiempo de trabajo del vendedor (convertido a costo por hora). Al final, se incluye una visualización del recorrido óptimo sobre el mapa de Colombia y un GIF que muestra paso a paso cómo se construye esa ruta.

Instalación y carga de paquetes

# Cargar paquetes necesarios
library(dplyr)
library(ggplot2)
library(gganimate)
library(leaflet)

Ciudades principales de Colombia

Creamos un data.frame con las 13 ciudades principales de Colombia y sus respectivas coordenadas geográficas (latitud y longitud).

ciudades <- data.frame(
  ciudad = c("Bogotá", "Medellín", "Cali", "Barranquilla", "Cartagena", 
             "Bucaramanga", "Pereira", "Manizales", "Ibagué", "Pasto",
             "Cúcuta", "Montería", "Santa Marta"),
  lat = c(4.7110, 6.2442, 3.4516, 10.9639, 10.3910,
          7.1193, 4.8133, 5.0703, 4.4389, 1.2136,
          7.8941, 8.7479, 11.2408),
  lon = c(-74.0721, -75.5812, -76.5319, -74.7964, -75.4794,
          -73.1227, -75.6961, -75.5138, -75.2322, -77.2811,
          -72.5050, -75.8814, -74.1990)
)

ciudades_nombres <- ciudades$ciudad
ciudades_cantidad <- nrow(ciudades)

Para obtener las matrices de distancias, de horas y de precios de peajes entre cada par de ciudades, se hace uso del siguiente servicio web que pertenece al Ministerio de Transporte de Colombia. Esta es una fuente muy confiable ya que es gubernamental.

Obtener la información de la anterior fuente fué la mejor alternativa debido a la confianza de los datos. Lamentablemente, no es una API donde se pueda obtener esta información a partir de consultas, tampoco se puede descargar la matriz con todos los datos. Para obtener estos datos, manualmente se tuvo que ingresar cada par de ciudades para obtener la información de la distancia, del tiempo promedio en el mejor de los casos y del costo de peajes.

Matriz de horas entre cada par de ciudades

horas_vector <- c(
  0, 5.83, 5.85, 12.03, 12.52, 5.83, 4.08, 4.2, 2.95, 9.42, 7.43, 9.95, 11.28, 5.83, 0, 5.20, 8.33, 7.53, 5.17, 2.63, 2.55, 4.12, 9.35, 7.85, 4.95, 9.40, 5.77, 5.20, 0, 13.55, 12.73, 9.03, 2.75, 3.27, 3.23, 4.57, 11.73, 10.17, 13.57, 12.03, 8.33, 13.55, 0, 1.6, 6.88, 10.98, 10.90, 11.52, 17.70, 8.30, 4.28, 1.33, 12.52, 7.53, 12.73, 1.60, 0, 7.38, 10.18, 10.10, 11.65, 16.88, 8.78, 3.15, 2.87, 5.83, 5.17, 9.03, 6.88, 7.38, 0, 6.48, 6.03, 6.23, 13.18, 2.68, 8.32, 6.13, 4, 2.63, 2.75, 10.98, 10.18, 6.48, 0, 0.7, 1.47, 6.09, 9.17, 7.60, 11.02, 4.20, 2.55, 3.27, 10.90, 10.10, 6.03, 0.70, 0, 2.18, 7.42, 8.73, 7.52, 10.57, 2.95, 4.20, 3.32, 11.52, 11.75, 6.23, 1.55, 2.27, 0, 7.47, 8.92, 9.17, 10.77, 9.42, 9.35, 4.93, 17.70, 16.88, 13.18, 6.90, 7.42, 7.38, 0, 15.47, 14.32, 17.72, 7.43, 7.85, 11.73, 8.30, 8.78, 2.68, 9.17, 8.73, 8.92, 15.47, 0, 9.72, 7.55, 9.95, 4.95, 10.17, 4.28, 3.15, 8.32, 7.60, 7.52, 9.08, 14.32, 9.72, 0, 5.33, 11.28, 9.40, 13.57, 1.33, 2.87, 6.13, 11.02, 10.57, 10.77, 17.72, 7.55, 5.33, 0
)

horas_matriz <- matrix(
  horas_vector, 
  nrow = ciudades_cantidad, 
  ncol = ciudades_cantidad, 
  byrow = TRUE,
  dimnames = list(ciudades_nombres, ciudades_nombres)
)

horas_matriz <- horas_matriz * 1.35
  • Se hace un aumento del 35 % en las horas que tomaría cada recorrido entre dos ciudades para tener en cuenta posibles retrazos por algún tipo de imprevisto o eventualidad.

Visualización parcial de la matriz

A continuación, mostramos una parte de la matriz de horas para validar los resultados.

round(horas_matriz[1:5, 1:5], 1)
##              Bogotá Medellín Cali Barranquilla Cartagena
## Bogotá          0.0      7.9  7.9         16.2      16.9
## Medellín        7.9      0.0  7.0         11.2      10.2
## Cali            7.8      7.0  0.0         18.3      17.2
## Barranquilla   16.2     11.2 18.3          0.0       2.2
## Cartagena      16.9     10.2 17.2          2.2       0.0

Matriz de distancias entre cada par de ciudades

distancias_viales_vector <- c(
  0.0, 417.42, 447.76, 971.37, 1038.66, 402.10, 309.80, 292.91, 196.75, 770.34, 577.67, 804.78, 927.16, 417.42, 0.0, 414.42, 695.13, 647.92, 398.48, 210.94, 196.73, 323.95, 793.31, 595.77, 402.49, 787.97, 439.94, 414.42, 0.0, 1109.55, 1062.34, 766.71, 207.98, 254.40, 250.98, 371.34, 964.00, 816.91, 1197.26, 971.37, 695.13, 1109.55, 0.0, 128.49, 575.06, 906.07, 891.87, 998.29, 1488.44, 71.23, 348.61, 98, 1038.66, 647.92, 1062.34, 128.49, 0.0, 642.35, 858.86, 844.66, 971.88, 1441.23, 777.52, 272.17, 230.72, 402.10, 398.48, 761.23, 575.06, 642.35, 0.0, 557.75, 511.80, 523.53, 1140.12, 197.29, 714.89, 530.85, 301.98, 210.94, 207.98, 906.07, 858.86, 557.75, 0.0, 50.3, 113.02, 586.87, 755.05, 613.43, 988.30, 292.91, 196.73, 254.40, 891.87, 844.66, 511.80, 50.93, 0.0, 163.94, 633.30, 709.09, 599.23, 942.35, 196.75, 331.77, 258.80, 998.29, 979.70, 523.53, 120.84, 171.76, 0.0, 637.69, 720.82, 734.27, 954.07, 770.34, 793.31, 384.17, 1488.44, 1441.23, 1145.60, 586.87, 633.30, 629.87, 0.0, 1347.89, 1195.80, 1576.15, 577.67, 595.77, 958.52, 710.23, 777.52, 197.29, 755.05, 709.09, 720.82, 1347.89, 0.0, 850.06, 666.02, 804.78, 402.49, 816.91, 348.61, 272.17, 714.89, 613.43, 599.23, 726.45, 1195.80, 850.06, 0.0, 441.44, 927.16, 787.97, 1191.78, 98, 230.72, 530.85, 988.30, 942.35, 954.07, 1570.67, 666.02, 441.44, 0.0
)

distancias_viales_matriz <- matrix(
  distancias_viales_vector, 
  nrow = ciudades_cantidad, 
  ncol = ciudades_cantidad, 
  byrow = TRUE,
  dimnames = list(ciudades_nombres, ciudades_nombres)
)

Visualización parcial de la matriz

A continuación, mostramos una parte de la matriz de distancias viales (en km) para validar los resultados.

round(distancias_viales_matriz[1:5, 1:5], 1)
##              Bogotá Medellín   Cali Barranquilla Cartagena
## Bogotá          0.0    417.4  447.8        971.4    1038.7
## Medellín      417.4      0.0  414.4        695.1     647.9
## Cali          439.9    414.4    0.0       1109.6    1062.3
## Barranquilla  971.4    695.1 1109.6          0.0     128.5
## Cartagena    1038.7    647.9 1062.3        128.5       0.0

Matriz de precios de peajes entre cada par de ciudades

peajes_vector <- c(
  0, 98400, 119600, 123800, 122300, 51500, 82900, 59700, 53100, 93600, 41200, 199600, 103600, 98400, 0, 138900, 145600, 133700, 68800, 82400, 66300, 112200, 151100, 90000, 101200, 165800, 119600, 138900, 0, 284500, 272600, 128100, 56500, 72600, 66500, 34300, 149300, 240100, 190400, 123800, 145600, 284500, 0, 23100, 94000, 228000, 211900, 144100, 296700, 76300, 57300, 20200, 122300, 133700, 272600, 23100, 0, 92500, 216100, 200000, 245900, 284800, 74800, 75500, 43300, 51500, 68800, 120500, 94000, 92500, 0, 64000, 47900, 61600, 132700, 21200, 105300, 73800, 82900, 82400, 56500, 228000, 216100, 64000, 0, 16100, 29800, 68700, 85200, 183600, 126300, 59700, 66300, 72600, 211900, 200000, 47900, 16100, 0, 45900, 84800, 69100, 167500, 110200, 53100, 112200, 66500, 144100, 245900, 61600, 29800, 45900, 0, 78700, 82800, 213400, 123900, 93600, 151100, 34300, 296700, 284800, 140300, 68700, 84800, 78700, 0, 146000, 252300, 202600, 41200, 90000, 141700, 76300, 74800, 21200, 85200, 69100, 82800, 146000, 0, 87600, 56100, 199600, 101200, 240100, 57300, 75500, 105300, 183600, 167500, 213400, 252300, 87600, 0, 77500, 103600, 165800, 182800, 20200, 43300, 73800, 126300, 110200, 123900, 195000, 56100, 77500, 0
)

peajes_matriz <- matrix(
  peajes_vector, 
  nrow = ciudades_cantidad, 
  ncol = ciudades_cantidad, 
  byrow = TRUE,
  dimnames = list(ciudades_nombres, ciudades_nombres)
)

Visualización parcial de la matriz

A continuación, mostramos una parte de la matriz de precios de peajes (en COP) para validar los resultados.

round(peajes_matriz[1:5, 1:5], 1)
##              Bogotá Medellín   Cali Barranquilla Cartagena
## Bogotá            0    98400 119600       123800    122300
## Medellín      98400        0 138900       145600    133700
## Cali         119600   138900      0       284500    272600
## Barranquilla 123800   145600 284500            0     23100
## Cartagena    122300   133700 272600        23100         0

Estimación de costos de viaje entre ciudades

Cálculo del costo por kilómetro

El costo por kilómetro se calcula con base en los siguientes supuestos:

  • Se asume que el recorrido se realiza en un Chevrolet Camaro 3.6L V6 modelo 2022.
  • Este vehículo tiene un consumo combinado estimado de 10.7 litros por cada 100 km, lo cual corresponde a un rendimiento de aproximadamente 9.35 km por litro.
  • Se considera un precio promedio de la gasolina de 16,300 COP por galón, lo que equivale a aproximadamente 4,308 COP por litro
    (1 galón = 3.785 litros).

Cálculo del costo por km:

\[ \text{Costo por km} = \frac{4,308\ \text{COP}}{9.35\ \text{km/L}} \approx \mathbf{460\ COP\ por\ km} \]

Este valor se utiliza como base para estimar el gasto en combustible entre cada par de ciudades.

Valor por hora del conductor

  • 60.000 COP

    • Se asume que trabaja 200 horas al mes y que esas 200 horas le valen 12.000.000 COP

  • Combustible = distancia * costo_km

costo_km <- 460           # pesos colombianos por kilómetro recorrido
valor_hora <- 60000       # valor hora del vendedor

Cálculo del costo total entre ciudades

Creamos una matriz con el costo total entre cada par de ciudades, donde:

\[ \text{costo_total} = (\text{tiempo} \cdot \text{valor_hora}) + \text{peajes} + (\text{distancia} \cdot \text{costo_km}) \]

En la diagonal principal, el costo será cero porque no hay desplazamiento.

costos_matriz <- matrix(0, nrow = ciudades_cantidad, ncol = ciudades_cantidad)
colnames(costos_matriz) <- ciudades_nombres
rownames(costos_matriz) <- ciudades_nombres

for (i in 1:ciudades_cantidad) {
  for (j in 1:ciudades_cantidad) {
    if (i != j) {
      distancia_vial <- distancias_viales_matriz[i, j]
      horas <- horas_matriz[i, j]
      costo_peaje <- peajes_matriz[i, j]
      combustible <- distancia_vial * costo_km
      
      costos_matriz[i, j] <- horas * valor_hora + costo_peaje + combustible
    }
  }
}

Visualización parcial de la matriz de costos

round(costos_matriz[1:5, 1:5], 0)
##               Bogotá Medellín    Cali Barranquilla Cartagena
## Bogotá             0   762643  799420      1545060   1614204
## Medellín      762643        0  750733      1140090   1041673
## Cali          789342   750733       0      1892443   1792406
## Barranquilla 1545060  1140090 1892443            0    211805
## Cartagena    1614204  1041673 1792406       211805         0

Observaciones

  • Esta matriz de costos será usada en los algoritmos de optimización: algoritmos genéticos y colonias de hormigas para minimizar el costo total del recorrido.

Algoritmo Genético para el Problema del Viajero (TSP)

Implementaremos un algoritmo genético (GA) para minimizar el costo total del recorrido entre ciudades.

Parámetros del algoritmo

generaciones_cantidad <- 200
poblacion_cantidad <- 100
tasa_mutacion <- 0.1

Función para evaluar una ruta

Dado un vector con el orden de las ciudades, esta función suma los costos entre ciudades consecutivas.

evaluar_ruta <- function(ruta, matriz_costos) {

  total <- 0

  for (i in 1:(ciudades_cantidad - 1)) {
    total <- total + matriz_costos[ruta[i], ruta[i + 1]]
  }

  # Regresar a la ciudad inicial
  total <- total + matriz_costos[ruta[ciudades_cantidad], ruta[1]]

  return(total)
}

Función para selección por torneo

seleccionar_padres_por_torneo <- function(poblacion, puntajes) {

  padres <- vector("list", poblacion_cantidad)
  
  for (i in 1:poblacion_cantidad) {
    torneo <- sample(1:poblacion_cantidad, 2)
    if (puntajes[torneo[1]] < puntajes[torneo[2]]) {
      padres[[i]] <- poblacion[[torneo[1]]]
    } else {
      padres[[i]] <- poblacion[[torneo[2]]]
    }
  }
  
  return(padres)
}

Funciones de cruzamiento y mutación

cruzar <- function(padre1, padre2) {
  punto <- sample(2:(ciudades_cantidad - 1), 1)
  hijo <- padre1[1:punto]
  hijo <- c(hijo, padre2[!padre2 %in% hijo])
  return(hijo)
}
mutar <- function(ruta) {
  if (runif(1) < tasa_mutacion) {
    pos <- sample(1:ciudades_cantidad, 2)
    ruta[pos] <- ruta[rev(pos)]
  }
  return(ruta)
}

Función para generar nueva población

generar_nueva_poblacion <- function(padres) {
  nueva_poblacion <- list()
  
  for (i in seq(1, poblacion_cantidad, by = 2)) {
    hijo1 <- mutar(cruzar(padres[[i]], padres[[i+1]]))
    hijo2 <- mutar(cruzar(padres[[i+1]], padres[[i]]))
    nueva_poblacion <- c(nueva_poblacion, list(hijo1, hijo2))
  }
  
  return(nueva_poblacion)
}

Implementación del algoritmo genético

# Inicializar población
poblacion <- replicate(poblacion_cantidad, sample(ciudades_nombres), simplify = FALSE)

for (gen in 1:generaciones_cantidad) {
  # Evaluar población
  puntajes <- sapply(poblacion, function(ruta) evaluar_ruta(ruta, costos_matriz))

  padres <- seleccionar_padres_por_torneo(poblacion, puntajes)
  poblacion <- generar_nueva_poblacion(padres)
}

Solución final

puntajes_finales <- sapply(poblacion, function(ruta) evaluar_ruta(ruta, costos_matriz))
mejor_ruta <- poblacion[[which.min(puntajes_finales)]]
mejor_costo <- min(puntajes_finales)

ruta_legible <- paste(c(mejor_ruta, mejor_ruta[1]), collapse = " → ")
cat("Mejor ruta encontrada:\n", ruta_legible, "\n\n")
## Mejor ruta encontrada:
##  Bogotá → Ibagué → Cali → Pasto → Pereira → Manizales → Medellín → Montería → Cartagena → Santa Marta → Barranquilla → Cúcuta → Bucaramanga → Bogotá
cat("Costo total: ", format(mejor_costo, big.mark = ","), "COP\n")
## Costo total:  6,286,303 COP

Visualización del recorrido en un mapa

Graficaremos el recorrido óptimo calculado por el Algoritmo Genetico usando coordenadas geográficas y líneas conectadas.

Preparación de datos

# Se obtienen los índices de las ciudades en el orden de mejor_ruta
indices <- match(mejor_ruta, ciudades_nombres)

# Se extraen y se reordenan las filas según mejor_ruta
ruta_coords <- ciudades[indices, ]

# Se agrega el punto inicial al final para cerrar el ciclo
ruta_coords <- rbind(ruta_coords, ruta_coords[1, ])

Mapa del recorrido óptimo

leaflet_output <- addCircleMarkers(
  addPolylines(
    addTiles(
      leaflet(ruta_coords)
    ),
    lng = ~lon, lat = ~lat,
    color = "blue", weight = 2,
    label = paste0("De ", ruta_coords$ciudad[-nrow(ruta_coords)],
                   " a ", ruta_coords$ciudad[-1])
  ),
  lng = ~lon, lat = ~lat,
  label = ~ciudad,
  color = "red", radius = 5,
  fillOpacity = 0.8
)

leaflet_output

Algoritmo de Colonias de Hormigas (ACO) para el TSP

Este algoritmo simula cómo las hormigas encuentran caminos cortos mediante feromonas y visibilidad (inverso de la distancia/costo).

###️ Parámetros del algoritmo

num_hormigas <- 50
num_iteraciones <- 100
alpha <- 1       # influencia de feromonas
beta <- 5        # influencia de visibilidad (1 / costo)
rho <- 0.5       # tasa de evaporación
Q <- 100         # cantidad de feromona depositada

Inicialización de estructuras

# Inicializar feromonas
feromonas <- matrix(1, nrow = ciudades_cantidad, ncol = ciudades_cantidad)
rownames(feromonas) <- ciudades_nombres
colnames(feromonas) <- ciudades_nombres

# Visibilidad (1 / costo), evitando división por cero
visibilidad <- 1 / (costos_matriz + diag(Inf, ciudades_cantidad))

Ciclo principal de ACO

mejor_ruta_aco <- NULL
mejor_costo_aco <- Inf

for (iter in 1:num_iteraciones) {
  rutas <- list()
  costos <- numeric(num_hormigas)
  
  for (k in 1:num_hormigas) {
    visitadas <- sample(ciudades_nombres, 1)
    while (length(visitadas) < ciudades_cantidad) {
      actual <- tail(visitadas, 1)
      no_visitadas <- setdiff(ciudades_nombres, visitadas)
      
      probabilidades <- sapply(no_visitadas, function(ciudad) {
        feromona <- feromonas[actual, ciudad]^alpha
        visib <- visibilidad[actual, ciudad]^beta
        feromona * visib
      })
      
      probabilidades <- probabilidades / sum(probabilidades)
      siguiente <- sample(no_visitadas, 1, prob = probabilidades)
      visitadas <- c(visitadas, siguiente)
    }
    rutas[[k]] <- visitadas
    costos[k] <- evaluar_ruta(visitadas, costos_matriz)
    
    if (costos[k] < mejor_costo_aco) {
      mejor_costo_aco <- costos[k]
      mejor_ruta_aco <- visitadas
    }
  }
  
  # Evaporación de feromonas
  feromonas <- (1 - rho) * feromonas
  
  # Depositar nuevas feromonas
  for (k in 1:num_hormigas) {
    ruta <- rutas[[k]]
    for (i in 1:(ciudades_cantidad - 1)) {
      a <- ruta[i]
      b <- ruta[i + 1]
      feromonas[a, b] <- feromonas[a, b] + Q / costos[k]
      feromonas[b, a] <- feromonas[b, a] + Q / costos[k]
    }
    # Conectar final con inicio
    feromonas[ruta[ciudades_cantidad], ruta[1]] <- feromonas[ruta[ciudades_cantidad], ruta[1]] + Q / costos[k]
    feromonas[ruta[1], ruta[ciudades_cantidad]] <- feromonas[ruta[1], ruta[ciudades_cantidad]] + Q / costos[k]
  }
}

Resultado del ACO

mejor_ruta_aco
##  [1] "Pasto"        "Cali"         "Ibagué"       "Bogotá"       "Bucaramanga" 
##  [6] "Cúcuta"       "Santa Marta"  "Barranquilla" "Cartagena"    "Montería"    
## [11] "Medellín"     "Manizales"    "Pereira"
mejor_costo_aco
## [1] 6267715

Creación del GIF animado

1. Preparación de los datos del recorrido

# Nos aseguramos de tener el mejor recorrido
ruta <- mejor_ruta_aco
coords <- ciudades[match(ruta, ciudades$ciudad), ]

# Agregar punto de retorno al inicio
coords <- rbind(coords, coords[1, ])

# Crear data.frame animable paso a paso
anim_data <- data.frame()
for (i in 2:nrow(coords)) {
  paso <- coords[1:i, ]
  paso$frame <- i - 1
  paso$orden <- 1:nrow(paso)
  anim_data <- rbind(anim_data, paso)
}

2. Creación de la animación punto por punto

p <- ggplot(anim_data, aes(x = lon, y = lat)) +
  borders("world", regions = "Colombia", fill = "gray95", colour = "gray60") +
  geom_path(aes(group = frame), color = "blue", linewidth = 1) +
  geom_point(color = "red", size = 3) +
  geom_text(aes(label = ciudad), vjust = -1, size = 3) +
  labs(title = "Recorrido óptimo - Paso: {frame}") +
  coord_fixed(1.3) +
  theme_minimal()

anim <- p + transition_manual(frames = frame)

anim_save(
  "recorrido_paso_a_paso.gif", 
  animation = anim, 
  renderer = gifski_renderer(),
  width = 832,   # ancho en píxeles
  height = 1080,  # alto en píxeles
  units = "px"
)
## `nframes` and `fps` adjusted to match transition

Referencias

Surjanovic, S., & Bingham, D. (s.f.). Rosenbrock Function. Simon Fraser University. https://www.sfu.ca/~ssurjano/rosen.html

Surjanovic, S., & Bingham, D. (s.f.). Schwefel Function. Simon Fraser University. https://www.sfu.ca/~ssurjano/schwef.html