REPORTE TRABAJO 1: SOLUCIÓN DE PROBLEMAS DE OPTIMIZACIÓN CON MÉTODOS HEURÍSTICOS

REDES NEURONALES Y ALGORITMOS BIOINSPIRADOS

Presentado por:

Leonardo Federico Corona Torres
David Escobar Ruiz
Johan Sebastian Robles Rincón
Sebastián Soto Arcila

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

Introducción

El presente trabajo trata algunos de los métodos más relevantes de optimización numérica basada en manipulaciones del gradiente de una función objetivo (función a optimizar) y también algoritmos heurísticos basados en el comportamiento de la naturaleza, como lo son los algoritmos evolutivos. Todos estos métodos pueden llegar a ser útiles en distintos contextos en los que se busca optimizar una función para algún fin, como encontrar la mejor ruta entre ciudades o optimizar una función de pérdida en un contexto de aprendizaje de máquina. Cada método tiene particularidades interesantes que se van a explorar de manera general por medio de una breve implementación y ejecución con dos funciones muy utiilizadas para probar algoritmos de optimización: la función de Rosenbrock y la función de Rastrigin.

Metodología

Para implementar y documentar la optimización se emplea R Markdown, combinando código R y texto explicativo. Se utilizan librerías para renderizar gráficos y animaciones, librerías para implementar los métodos heurísticos y librerías para manipular los datos.

# DESCOMENTAR PARA INSTALAR PAQUETES FALTANTES
#install.packages("GA")
#install.packages("tidyverse")
#install.packages("viridis")
#install.packages("ggplot2")
#install.packages("gganimate")
#install.packages("gifski")
#install.packages("dplyr")
#install.packages("tidyr")
#install.packages("purrr")
#install.packages("pso")
#install.packages("leaflet")
#install.packages("osrm")
#install.packages("sf")
#install.packages("gor")

# DESCOMENTAR PARA INSTALAR PAQUETES FALTANTES
#install.packages("GA")
#install.packages("tidyverse")
#install.packages("viridis")
#install.packages("ggplot2")
#install.packages("gganimate")
#install.packages("gifski")
#install.packages("dplyr")
#install.packages("tidyr")
#install.packages("purrr")
#install.packages("pso")
#install.packages("leaflet")
#install.packages("osrm")
#install.packages("sf")
#install.packages("gor")

Para los casos 3D (tres variables) es difícil visualizar directamente la función de 3 dimensiones (espacio de 4 variables), por lo que se mostrarán únicamente los puntos óptimos alcanzados.

Parte 1. Optimización numérica

1.1 Selección e implementación de las funciones de prueba

Para observar el comportamiento de cada uno de los métodos de optimización implementados, se optó por la selección de dos funciones de prueba estándar para evaluar algoritmos de optimización: Función de Rosenbrock y Función de Rastrigin. Debido a su complejidad, ambas funciones son indicadores útiles para comparar la eficacia de optimizadores como los AG.

Previo a la implementación de los algoritmos en el lenguaje R y previo a la recopilación de resultados y la derivación de conclusiones a partir de estos, el equipo decidió realizar una breve investigación general de lo que son las funciones de prueba para problemas de optimización. La investigación se realizó de tal forma que, una vez concluída, se pudieran comprender los conceptos hasta el punto de estar en las capacidades responder las siguientes preguntas implícitamente en una breve sección de este trabajo:

• ¿Qué es una función de prueba en optimización?

• ¿Para qué se utiliza una función de prueba en optimización?

• ¿Cómo se pueden clasificar las funciones de prueba?

• ¿Cuáles son algunas de las funciones de prueba más utilizadas?

A continuación se presenta dicha sección para luego continuar con el estudio más detallado de las funciones escogidas.

1.1.1 Breve explicación de las funciones de prueba para problemas de optimización

Según (Yang, 2010), una función de prueba es una función con unas propiedades especiales que permite probar si el rendimiento de un método de optimización implementado es aceptable bajo las condiciones especiales que impone la función de prueba. 

Esto es especialmente útil para verificar que el método sea eficiente bajo distintas condiciones en las que se espera que se implemente, como por ejemplo, casos en los que la función tiene múltiples mínimos y/o máximos locales.

Según (Molga, 2005), las funciones de prueba se pueden ubicar en una de las siguientes clases, todas siendo funciones continuas:

• Clase 1: Unimodal, convexa, multidimensional.
• Clase 2: Multimodal, dos dimensiones con un número pequeño de extremos locales.
• Clase 3: Multimodal, dos dimensiones con un gran número de extremos locales.
• Clase 4: Multimodal, multidimensional, con un número amplio de extremos locales.

En el caso de las funciones elegidas, la función de Rosenbrock se clasificaría como Clase 3 y la de Rastrigin como Clase 2.

Como ejemplo, en la Figura 1 se presentan algunas funciones de prueba que no se eligieron para este trabajo, pero que son igual de relevantes, importantes y comúnmente utilizadas en la práctica (Molga, 2005):

• Función de De Jong (Clase 1).

• Función de Griewangk (Clase 2).

• Función de Langermann (Clase 3).

• Función de Ackley (Clase 4).

1.1.2 Implementación de funciones de graficación

Se definen las funciones a usar para generar los gráficos de cada función de esta sección:

generate_function_graph <- function(func_2d, low_bound, upper_bound) {
  x1 <- seq(low_bound, upper_bound, length.out = 50)
  x2 <- seq(low_bound, upper_bound, length.out = 50)
  f_x <- outer(x1,x2, FUN = func_2d)
  
  colores        <- viridis::magma(n = 100, alpha = 0.7)
  z.facet.center <- (f_x[-1, -1] + f_x[-1, -ncol(f_x)] +
                       f_x[-nrow(f_x), -1] +
                       f_x[-nrow(f_x), -ncol(f_x)])/4
  z.facet.range  <- cut(z.facet.center, 100)
  
  par(mai = c(0,0,0,0))
  persp(x = x1, y = x2, z = f_x,
        shade = 0.8,
        phi = 30,
        theta = 30,
        col = colores[z.facet.range],
        axes = FALSE)
}

generate_level_curves <- function(func_2d, low_bound, upper_bound) {
  n_length <- 100
  x1 <- seq(low_bound, upper_bound, length.out = n_length)
  x2 <- seq(low_bound, upper_bound, length.out = n_length)
  X <- expand.grid(x1, x2)
  z <- func_2d(X[,1], X[,2])
  Z <- matrix(z, ncol = n_length, nrow = n_length)
  contour(
    x = x1,
    y = x2,
    z = Z,
    nlevels = 100,
    las = 1,
    xlab = expression(x[1]),
    ylab = expression(x[2]),
    main = expression(paste(
      "Función de Rosenbrock: ",
      f(x[1],x[2])==100*(x[2]-x[1]^2)^2+(1-x[1])^2
    )),
    sub = "Curvas de nivel de la función"
  )
}

generate_level_curves_2 <- function(func_2d, low_bound, upper_bound) {
  x1 <- seq(low_bound, upper_bound, length.out = 50)
  x2 <- seq(low_bound, upper_bound, length.out = 50)
  
  datos <- expand.grid(x1 = x1, x2 = x2)
  datos <- datos %>%
           mutate(f_x = map2_dbl(x1, x2, .f = func_2d))
  
  f_level_curves <- ggplot(data = datos, aes(x = x1, y = x2, z = f_x)) +
    geom_contour(aes(colour = stat(level)), bins = 30) +
    labs(title = "f(x1,x2) = x1^2 + x2^2") +
    theme_bw() +
    theme(legend.position = "none")
  return(f_level_curves)
}

generate_countour <- function(func_2d, low_bound, upper_bound) {
  n_length <- 100
  x1 <- seq(low_bound, upper_bound, length.out = n_length)
  x2 <- seq(low_bound, upper_bound, length.out = n_length)
  expand.grid(X1 = x1, X2 = x2) %>%
    mutate(Z = func_2d(X1, X2)) %>%
    ggplot(aes(X1, X2, z = Z)) +
    geom_contour() +
    geom_contour_filled()
}

1.1.3 Función de Rosenbrock

Es una función no convexa introducida por Rosenbrock en 1960​ (Wikipedia, s.f., Rosenbrock function). Su paisaje forma un valle curvo estrecho que dificulta la convergencia hacia el mínimo.

Definición en n dimensiones:

\(f(X) = \sum_{i=1}^{n-1}100(X_{i+1}-X_i^2)^2 + (1-X_i)^2\)

Definición en 2D:

\(f(x,y)= 100(y-x^2)^2 + (1-x)^2\)

Definición en 3D:

\(f(x,y,z)=100[(y-x^2)^2+(z-y^2)^2] + (1-x)^2 + (1-y)^2\)

Algunas características y propiedades importantes  son las siguientes:

• Mínimo global:

  \(x_1,...,x_n =1, f(x_1,...,x_n) = 0\)

• Dominio de búsqueda: 

  \(-\infty \leq X_i \leq \infty\)

  \(1 \leq i \leq n\)

• Conocida también como la función banana de Rosenbrock.

• Tiene forma de valle, el cual es trivial encontrarlo. Sin embargo, la convergencia al mínimo global es difícil.

Estas son algunas hipótesis y expectativas que se tuvieron para el rendimiento de cada método implementado con esta función:

• Método de descenso de gradiente:

  • Si se ubica el punto inicial sobre las rectas tangentes con mayor gradiente a la función entonces el método llegará más rápido a un mínimo.
  • Por el contrario, si la recta tangente tiene una gradiente más cercana a cero, el método llegará más lentamente, es decir, requerirá de más iteraciones para llegar al punto mínimo.

• Método de algoritmos evolutivos:

  • Los criterios para seleccionar los individuos a reproducir va a ser muy importante para el rendimiento del método.

• Método de optimización de partículas:

  • Este método no va a tener muchas dificultades para encontrar el punto óptimo es la función de Rosenbrock.
  • El criterio para parar las iteraciones no va a cambiar mucho los resultados en este caso.

• Método de evolución diferencial:

  • Se comporta similar a los métodos evolutivos en 2 dimensiones.

La implementación de la función de Rosenbrock en 2, 3 y N dimensiones se muestra a continuación:

generate_countour(f_rosenbrock_2d, -5, 5)

Las gráficas de la función de Rosenbrock en 2 y 3 dimensiones se muestran a continuación:

generate_countour(f_rosenbrock_2d, -5, 5)

generate_level_curves(f_rosenbrock_2d, -5, 5)

generate_function_graph(f_rosenbrock_2d, -5, 5)

1.1.4 Función de Rastrigin

La función Rastrigin (1974) es una función no convexa y altamente multimodal, con numerosos mínimos locales, lo que la hace difícil de optimizar​ (Wikipedia, s.f., Rastrigin function).

Presenta múltiples mínimos locales dispuestos de forma regular​, lo que lo convierte en un desafío típico para algoritmos de optimización.

Definición en n dimensiones:

\(f(X)=An + \sum_{i=1}^n[X_i^2 -A\cos(2\pi X_i)], A=10\)

Definición en 2D:

\(f(x,y)=x^2+y^2+A[2-\cos(2\pi x) - \cos(2\pi y)], A=10\)

Definición en 3D:

\(f(x,y)=x^2+y^2+z^2+A[3-\cos(2\pi x) - \cos(2\pi y)-\cos(2\pi z)], A=10\)

Algunas características y propiedades importantes son las siguientes:

• Mínimo global:

  • \(f(0,...0)=0\)

• Dominio de búsqueda:

  • \(-5.12 \leq X_i \leq 5.12\)

• Particularmente, hallar el mínimo de esta función es un problema difícil debido a la larga cantidad de mínimos locales.

Estas son algunas hipótesis y expectativas que se tuvieron para el rendimiento de cada método implementado con esta función:

• Método de descenso de gradiente:

  • El éxito del método en alcanzar un mínimo dependerá enormemente del punto inicial elegido debido a los múltiples mínimos locales que tiene esta función.

• Método de algoritmos evolutivos:

  • Es importante hacer una buena representación de los individuos, ya que el rendimiento puede variar dependiendo de la selección.
  • La población inicial debe ser representativa, diferentes entre ellos y que no sean demasiados. Para esta función en particular encontrar la solución óptima va a depender mucho de la probabilidad de mutación

• Método de optimización de partículas:

  • La ubicación de la nube de partículas va a afectar el rendimiento en esta función, ya que el algoritmo puede quedarse estancado en un mínimo local si la ubicación no es la mejor.
  • El criterio para detener el algoritmo debe estar bien ajustado, para evitar estancamientos.

• Método de evolución diferencial:

  • Es más eficiente que los otros métodos evolutivos ya que su fuerte son las múltiples dimensiones.

La implementación de la función de Rastrigin en 2, 3 y N dimensiones se muestra a continuación:

# 2 Dimensiones 
f_rastrigin_2d <- function(x, y) {
  A = 10   
  f_value <- x^2 + y^2 + A*(2 - cos(2*pi*x) - cos(2*pi*y))   
  return(f_value) 
}  

# 3 Dimensiones 
f_rastrigin_3d <- function(x, y, z) {   
  A = 10   
  f_value <- x^2 + y^2 + z^2 + A*(3 - cos(2*pi*x) - cos(2*pi*y) - cos(2*pi*z))   
  return(f_value) 
}

# N Dimensiones
f_rastrigin <-function(x){
  A <- 10
  n <- length(x)
  z <- (A*n) + sum(x^2 - A*cos(2*pi*x))
  return(z)
}

Las gráficas de la función de Rastrigin en 2 y 3 dimensiones se muestran a continuación:

generate_countour(f_rastrigin_2d, -10, 10)

generate_level_curves(f_rastrigin_2d, -10, 10)

generate_level_curves_2(f_rastrigin_2d, -10, 10)

generate_function_graph(f_rastrigin_2d, -10, 10)

1.3 Método de descenso de gradiente

1.3.1 Implementación en R de descenso por gradiente

Implementación de derivada parcial

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

Implementación del gradiente

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)
}

Implementación de derivada del gradiente

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

Implementación de matriz Hessiana

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

Implementación completa de optimizador multivariado por descenso de gradiente

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

1.3.2 Optimización de la función de Rosenbrock en 2 dimensiones

# Ejecución del método
sol_rosen <- optimizador_mult_numdev(f_rosenbrock,x0=c(-4,-4),eta=1)

# Graficación del proceso de optimización
n_length <- 100
x1 <- seq(-5, 5, length.out = n_length)
x2 <- seq(-5, 5, length.out = n_length)
X <- expand.grid(x1, x2)
z <- f_rosenbrock_2d(X[,1], X[,2])
Z <- matrix(z, ncol = n_length, nrow = n_length)
contour(
  x = x1,
  y = x2,
  z = Z,
  nlevels = 100,
  las = 1,
  xlab = expression(x[1]),
  ylab = expression(x[2]),
  main = expression(paste(
    "Función de Rosenbrock: ",
    f(x[1],x[2])==100*(x[2]-x[1]^2)^2+(1-x[1])^2)
  ),
  sub = "Curvas de nivel de la función"
)
lines(sol_rosen, type="b",cex=1.5,col="red")

Animación de la optimización

# Transforma la solución en un data.frame
sol <- as.data.frame(sol_ras2d)
colnames(sol) <- c("x1", "x2")
sol$step <- 1:nrow(sol)

# Prepara los datos de fondo
x1 <- seq(-5, 5, length.out = 100)
x2 <- seq(-5, 5, length.out = 100)
X <- expand.grid(x1 = x1, x2 = x2)

X$z <- f_rosenbrock_2d(X$x1, X$x2)

# Crea la animación
p <- ggplot() +
  geom_contour(data = X, aes(x = x1, y = x2, z = z), bins = 20, color = "gray70") +
  geom_point(data = sol, aes(x = x1, y = x2), color = "red", size = 3) +
  geom_path(data = sol, aes(x = x1, y = x2), color = "red", linewidth = 1) +
  transition_reveal(step) +
  labs(
    title = "Optimización sobre la función de Rosenbrock con descenso del gradiente",
    x = expression(x[1]),
    y = expression(x[2])
  ) +
  theme_minimal()

animate(p, fps = 10, duration = 5, width = 600, height = 500, renderer = gifski_renderer("rosenbrock_opt.gif"))

1.3.3 Optimización de la función de Rosenbrock en 3 dimensiones

# Ejecución del método
sol_rosen <- optimizador_mult_numdev(f_rosenbrock,x0=c(-4,-4,-4),eta=1)

# Graficación del proceso de optimización
n_length <- 100
x1 <- seq(-5, 5, length.out = n_length)
x2 <- seq(-5, 5, length.out = n_length)
x3 <- seq(-5, 5, length.out = n_length)
X <- expand.grid(x1, x2, x3)
z <- f_rosenbrock_3d(X[,1], X[,2], X[,3])
Z <- matrix(z, ncol = n_length, nrow = n_length)
contour(
  x = x1,
  y = x2,
  z = Z,
  nlevels = 100,
  las = 1,
  xlab = expression(x[1]),
  ylab = expression(x[2]),
  main = expression(paste(
    "Función de Rosenbrock: ",
    f(x[1],x[2])==100*(x[2]-x[1]^2)^2+(1-x[1])^2)
  ),
  sub = "Curvas de nivel de la función"
)
lines(sol_rosen, type="b",cex=1.5,col="red")

1.3.4 Optimización de la función de Rastrigin en 2 dimensiones

sol_ras2d <- optimizador_mult_numdev(f_rastrigin,x0=c(4.5,4.5),eta=2)

n_length <- 100
x1 <- seq(-5.12, 5.12, length.out = n_length)
x2 <- seq(-5.12, 5.12, length.out = n_length)
X <- expand.grid(x1, x2)
z <- f_rastrigin_2d(X[,1], X[,2])
Z <- matrix(z, ncol = n_length, nrow = n_length)
contour(
  x = x1,
  y = x2,
  z = Z,
  nlevels = 10,
  las = 1,
  xlab = expression(x[1]),
  ylab = expression(x[2]),
  main = expression(paste(
    "Función de Rastrigin: ",
    f(x[1],x[2])==20 + x[1]^2 - 10*cos(2*pi*x[1]) + x[2]^2 - 10*cos(2*pi*x[2])
  )),
  sub = "Curvas de nivel de la función"
)
lines(sol_ras, type="b",cex=1.5,col="red")

Animación de la optimización

# Transforma la solución en un data.frame
sol <- as.data.frame(sol_ras2d)
colnames(sol) <- c("x1", "x2")
sol$step <- 1:nrow(sol)

#Prepara los datos de fondo
x1 <- seq(-5.12, 5.12, length.out = 100)
x2 <- seq(-5.12, 5.12, length.out = 100)
X <- expand.grid(x1 = x1, x2 = x2)

X$z <- f_rastrigin_2d(X$x1, X$x2)

# Crea la animación
p <- ggplot() +
  geom_contour(data = X, aes(x = x1, y = x2, z = z), bins = 20, color = "gray70") +
  geom_point(data = sol, aes(x = x1, y = x2), color = "red", size = 3) +
  geom_path(data = sol, aes(x = x1, y = x2), color = "red", linewidth = 1) +
  transition_reveal(step) +
  labs(
    title = "Optimización sobre la función de Rastrigin con descenso del gradiente",
    x = expression(x[1]),
    y = expression(x[2])
  ) +
  theme_minimal()

animate(p, fps = 10, duration = 5, width = 600, height = 500, renderer = gifski_renderer("rastrigin_opt.gif"))

1.3.5 Optimización de la función de Rastrigin en 3 dimensiones

# Ejecución del método
sol_ras3d <- optimizador_mult_numdev(f_rastrigin,x0=c(-4,-4,-4),eta=3)

# Graficación del proceso de optimización
n_length <- 100
x1 <- seq(-5.12, 5.12, length.out = n_length)
x2 <- seq(-5.12, 5.12, length.out = n_length)
x3 <- seq(-5.12, 5.12, length.out = n_length)
X <- expand.grid(x1, x2, x3)
z <- f_rastrigin_3d(X[,1], X[,2], X[,3])
Z <- matrix(z, ncol = n_length, nrow = n_length)
contour(
  x = x1,
  y = x2,
  z = Z,
  nlevels = 10,
  las = 1,
  xlab = expression(x[1]),
  ylab = expression(x[2]),
  main = expression(paste(
    "Función de Rastrigin: ",
    f(x[1],x[2])==20 + x[1]^2 - 10*cos(2*pi*x[1]) + x[2]^2 - 10*cos(2*pi*x[2])
  )),
  sub = "Curvas de nivel de la función"
)
lines(sol_rosen, type="b",cex=1.5,col="red")

1.4 Método de evolución diferencial

La evolución diferencial es un algoritmo de optimización inspirado en la evolución biológica. Funciona manteniendo una población de soluciones, y mejorándolas generación tras generación mediante operaciones de mutación, recombinación y selección.

• Mutación: se combinan 3 individuos distintos de la población para crear una variante.
• Recombinación: se mezcla esa variante con el individuo actual.
• Selección: se escoge el mejor entre el original y el nuevo.

Este proceso se repite varias veces hasta encontrar una solución óptima.

1.4.1 Implementación en R de evolución diferencial

evolucion_diferencial <- function(fun_obj, dim = 2, NP = 30, F = 0.8, CR = 0.9,
                                  gens = 100, bounds = c(-5, 5)) {

  # Inicializar población
  poblacion <- matrix(runif(NP * dim, bounds[1], bounds[2]), ncol = dim)
  fitness <- apply(poblacion, 1, fun_obj)

  historial <- numeric(gens)
  mejores <- matrix(NA, gens, dim)

  for (gen in 1:gens) {
    for (i in 1:NP) {
      # Seleccionar 3 índices distintos
      indices <- sample(setdiff(1:NP, i), 3)
      x1 <- poblacion[indices[1], ]
      x2 <- poblacion[indices[2], ]
      x3 <- poblacion[indices[3], ]

      # Mutación
      mutado <- x1 + F * (x2 - x3)

      # Recombinar
      trial <- poblacion[i, ]
      jrand <- sample(1:dim, 1)
      for (j in 1:dim) {
        if (runif(1) < CR || j == jrand) {
          trial[j] <- mutado[j]
        }
      }

      # Selección
      if (fun_obj(trial) < fitness[i]) {
        poblacion[i, ] <- trial
        fitness[i] <- fun_obj(trial)
      }
    }

    # Guardar mejor resultado
    best_idx <- which.min(fitness)
    historial[gen] <- fitness[best_idx]
    mejores[gen, ] <- poblacion[best_idx, ]
  }

  list(mejor = poblacion[which.min(fitness), ],
       valor = min(fitness),
       historial = historial,
       trayectoria = mejores)
}

1.4.2 Optimización de la función de Rosenbrock en 2 dimensiones

# Ejecución del método
set.seed(123)
res_rosen <- evolucion_diferencial(f_rosenbrock, dim = 2)

# Mostrar mejor solución
res_rosen$mejor
res_rosen$valor

# Graficar trayectoria
x1 <- seq(-3, 3, length.out = 100)
x2 <- seq(-3, 3, length.out = 100)
z <- outer(x1, x2, Vectorize(function(x, y) f_rosenbrock(c(x, y))))
contour(x1, x2, z, nlevels = 50,
        main = "Rosenbrock 2D - Trayectoria", xlab = "x", ylab = "y")
lines(res_rosen$trayectoria[,1], res_rosen$trayectoria[,2], col = "red", type = "b")

En este gráfico se muestran las curvas de nivel de la función de Rosenbrock en dos dimensiones. Estas curvas representan líneas donde la función tiene igual valor, y el valle curvado al centro es donde está el mínimo global (en el punto (1,1)(1,1)).

La línea roja representa la trayectoria que siguió el algoritmo de evolución diferencial durante las iteraciones. Se puede observar cómo el enjambre de soluciones se fue acercando progresivamente hacia el mínimo, mejorando su posición en cada generación.

El resultado final obtenido fue: [1] 1.000000 1.000001 Valor mínimo encontrado: 1.91e-12

1.4.3 Optimización de la función de Rosenbrock en 3 dimensiones

# Ejecución del método
res_rosen_3d <- evolucion_diferencial(f_rosenbrock, dim = 3)
res_rosen_3d$mejor
res_rosen_3d$valor

En esta parte del informe se muestran los resultados numéricos para la optimización de las funciones en 3D. Dado que no se puede visualizar fácilmente en una gráfica 3D de trayectoria, se reportan las mejores posiciones y valores obtenido.

1.4.4 Optimización de la función de Rastrigin en 2 dimensiones

# Ejecución del método
set.seed(456)
res_ras <- evolucion_diferencial(f_rastrigin, dim = 2)

# Mostrar mejor solución
res_ras$mejor
res_ras$valor

# Graficar trayectoria
x1 <- seq(-5.12, 5.12, length.out = 100)
x2 <- seq(-5.12, 5.12, length.out = 100)
z <- outer(x1, x2, Vectorize(function(x, y) f_rastrigin(c(x, y))))
contour(x1, x2, z, nlevels = 50,
        main = "Rastrigin 2D - Trayectoria", xlab = "x", ylab = "y")
lines(res_ras$trayectoria[,1], res_ras$trayectoria[,2], col = "blue", type = "b")

Aquí se grafican las curvas de nivel de la función de Rastrigin, que es multimodal, es decir, tiene muchos mínimos locales (patrón ondulado). La búsqueda es mucho más compleja que en Rosenbrock.

La línea azul muestra cómo la evolución diferencial se mueve por el espacio de búsqueda y logra escapar de los mínimos locales hasta acercarse al óptimo global, que se encuentra en (0,0)(0,0).

Resultado obtenido: [1] 2.176697e-06 -2.015785e-07 Valor mínimo: 9.48e-10

1.4.5 Optimización de la función de Rastrigin en 3 dimensiones

# Ejecución del método
res_ras_3d <- evolucion_diferencial(f_rastrigin, dim = 3)
res_ras_3d$mejor
res_ras_3d$valor

En este caso, el algoritmo encontró una solución cercana al mínimo global, aunque no exacta. Esto es esperable, ya que Rastrigin es mucho más difícil en 3D debido a la gran cantidad de mínimos locales.

1.4.6 Conclusiones método de evolución diferencial

El algoritmo logró encontrar soluciones muy cercanas al mínimo global en las funciones de Rosenbrock y Rastrigin, tanto en 2D como en 3D.

En la función de Rosenbrock, se observó una convergencia estable y precisa, especialmente en dos dimensiones.

En la función de Rastrigin, que presenta muchos mínimos locales, el método evitó caer en estos y alcanzó buenos resultados.

El rendimiento se mantuvo sólido en tres dimensiones, aunque con una ligera pérdida de precisión en Rastrigin 3D debido a su complejidad.

Fue más robusto que el descenso por gradiente en funciones multimodales, ya que no necesita derivadas ni depende del punto inicial.

La evolución diferencial demostró ser un método confiable y efectivo para resolver problemas de optimización continua.

1.5 Método de optimización de partículas

El método de optimización por enjambre de partículas (PSO, por sus siglas en inglés) es un algoritmo inspirado en el comportamiento colectivo de animales como bandadas de aves o bancos de peces. Funciona mediante un conjunto de partículas (soluciones potenciales) que exploran el espacio de búsqueda moviéndose en función de su propia experiencia y la de sus vecinas. Cada partícula ajusta su posición y velocidad iterativamente para acercarse a la mejor solución conocida, guiada por su mejor posición histórica y la mejor posición global encontrada por el enjambre. Con el tiempo, las partículas tienden a converger hacia una solución óptima o cercana al óptimo.

1.5.1 Implementación en R de optimización de partículas

Se utiliza el paquete pos para implementar el método de la optimización de partículas. Adicionalmente, se implementa una función adicional para crear las animaciones del método de optimización.

list_to_matrix <- function(data) {
  for (i in 1:length(data)) {
    data[[i]] <- matrix(data[[i]], nrow=2, ncol=12)
  }
  return(data)
}

particle_swarm_optimization <- function(n,func,lower_bounds,upper_bounds) {
  set.seed(2001)
  o_min <- psoptim(rep(NA,n), func, lower=lower_bounds,upper=upper_bounds,control=list(fnscale=1e-8,trace=1,trace.stats=TRUE))
  o_max <- psoptim(rep(NA,n), func, lower=lower_bounds,upper=upper_bounds,control=list(fnscale=-1*(1e-8),trace=1,trace.stats=TRUE,s=30))
  print("=====================SUMMARY=====================")
  print("MINIMIZATION:")
  print("Point:")
  show(o_min$par)
  print("Value:")
  show(func(o_min$par))
  print("MAXIMIZATION")
  print("Point:")
  show(o_max$par)
  print("Value:")
  show(func(o_max$par))
  return(list("o_min"=o_min, "o_max"=o_max))
}

animate_pso <- function(particles_positions, func, gif_name) {
  positions_per_iteration <- list_to_matrix(particles_positions)
  df <- map2_dfr(
    positions_per_iteration,
    .y = seq_along(positions_per_iteration),
    .f = function(mat, iter) {
      tibble(
        particle_id = 1:ncol(mat),
        x = mat[1, ],
        y = mat[2, ],
        iter = iter
      )
    }
  )
  x_seq <- seq(-10, 10, length.out = 100)
  y_seq <- seq(-10, 10, length.out = 100)
  grid <- expand.grid(x = x_seq, y = y_seq)
  grid$z <- with(grid, func(x, y))
  
  p <- ggplot() +
    geom_contour(data = grid, aes(x=x, y=y, z=z),
                 bins = 30, color = "gray") +
    geom_point(data=df, aes(x = x, y = y), color="red", size=2) +
    xlim(-10, 10) + ylim(-10, 10) +  # Adjust limits to your data range
    theme_minimal() +
    transition_manual(frames = iter) +
    labs(title = "Iteration: {current_frame}")
  #p + transition_reveal(agno)
  animate(p, fps = 5, renderer=gifski_renderer())
  anim_save(gif_name,p)
}

1.5.2 Optimización de la función de Rosenbrock en 2 dimensiones

# Parámetros a utilizar
n <- 2
lower_bounds <- -20
upper_bounds <- 20

# Ejecución del método
o_pso <- particle_swarm_optimization(n,f_rosenbrock,lower_bounds,upper_bounds)

# Graficación del proceso de optimización (minimización)
o_min <- o_pso$o_min
o_min_particles_positions <- o_min$stats$x
gif_name <- "pso_rosenbrock_min.gif"
animate_pso(o_min_particles_positions, f_rosenbrock_2d, gif_name)

# Graficación del proceso de optimización (maximización)
o_max <- o_pso$o_max
o_max_particles_positions <- o_max$stats$x
gif_name <- "pso_rosenbrock_max.gif"
animate_pso(o_max_particles_positions, f_rosenbrock_2d, gif_name)

1.5.3 Optimización de la función de Rosenbrock en 3 dimensiones

# Parámetros a utilizar
n <- 3
lower_bounds <- -5
upped_bounds <- 5

# Ejecución del método
o_pso <- particle_swarm_optimization(n,f_rosenbrock,lower_bounds,upper_bounds)
# Graficación del proceso de optimización

1.5.4 Optimización de la función de Rastrigin en 2 dimensiones

# Parámetros a utilizar
n <- 2
lower_bounds <- -5
upped_bounds <- 5

# Ejecución del método
o_pso <- particle_swarm_optimization(n,f_rastrigin,lower_bounds,upper_bounds)

# Graficación del proceso de optimización (minimización)
o_min <- o_pso$o_min
o_min_particles_positions <- o_min$stats$x
gif_name <- "pso_rastrigin_min.gif"
animate_pso(o_min_particles_positions, f_rastrigin_2d, gif_name)

# Graficación del proceso de optimización (maximización)
o_max <- o_pso$o_max
o_max_particles_positions <- o_max$stats$x
gif_name <- "pso_rastrigin_max.gif"
animate_pso(o_max_particles_positions, f_rastrigin_2d, gif_name)

1.5.5 Optimización de la función de Rastrigin en 3 dimensiones

# Parámetros a utilizar
n <- 3
lower_bounds <- -5
upped_bounds <- 5
# Ejecución del método
o_pso <- particle_swarm_optimization(n,f_rastrigin,lower_bounds,upper_bounds)

1.5.6 Conclusiones método de optimización de partículas

En el caso de la optimización maximizando las funciones, llegan muy rápido a una solución óptima (dentro de 10 iteraciones). Por otro lado, la optimización minimizando las funciones tardaba un poco más en llegar al óptimo, pero se acercaba mucho al mínimo global.

Estos resultados pueden deberse a la distribución de las partículas por el espacio de búsqueda, permitiendo una optimización más “abierta” a comparación con el método de descenso de gradiente.

1.6 Método de algoritmos evolutivos

Los algoritmos genéticos (AG) son técnicas de búsqueda heurística basadas en procesos de evolución natural​ . En un AG típico se define una función FITNESS que evalúa la calidad de cada solución candidata (individuo). A partir de una población inicial aleatoria, se iteran ciclos donde se seleccionan individuos más aptos, se combinan sus “genes” mediante cruces (crossover) y se introducen modificaciones aleatorias (mutaciones). Estos operadores evolucionan la población hacia regiones con mejor fitness.

Según (Scrucca, 2013), los GAs han sido exitosos en optimizar funciones continuas (diferenciables o no) y discretas. Entre los operadores genéticos clave se destacan:

• Selección: elige individuos con mayor fitness para reproducirse, imitando la supervivencia del más apto.

• Cruce (crossover): combina partes de dos soluciones parentales para generar descendencia, explorando nuevas regiones del espacio de búsqueda.

• Mutación: altera aleatoriamente parte de un individuo (por ejemplo, cambiando un valor de su vector de variables) para introducir diversidad genética y evitar estancamiento en óptimos locales.

Los algoritmos genéticos (AG) son metaheurísticas inspiradas en procesos evolutivos biológicos, que han demostrado eficacia en la búsqueda global de óptimos en funciones complejas​jstatsoft.org Los AG simulan la selección natural, la recombinación (cruce) y la mutación para iterativamente mejorar un conjunto de soluciones candidatas (población). Estas técnicas estocásticas son adecuadas para funciones no lineales, discontinuas o con múltiples óptimos locales donde los métodos basados en derivadas pueden fallar. Para evaluar la robustez de los GA, se realizarán múltiples ejecuciones independientes y se analizará la dispersión del fitness resultante.

Adicionalmente, se definen algunas funciones para poder mostrar por medio de una animación el proceso de optimización para las funciones de Rosenbrock y de Rastrigin.

animate_ga_optimization <- function(func) {
  # 1. Crear el entorno para almacenar la evolución de la población
  pop_data <- data.frame()
  
  # 2. Ejecutar el algoritmo genético, capturando las poblaciones
  ga_rastrigin <- ga(
    type = "real-valued",
    fitness = function(x) -func(x),
    lower = c(-5.12, -5.12), upper = c(5.12, 5.12),
    popSize = 50, maxiter = 50, run = 50,
    monitor = function(obj) {
      gen <- obj@iter
      pop <- obj@population
      df <- data.frame(
        X1 = pop[, 1],
        X2 = pop[, 2],
        Generacion = gen
      )
      pop_data <<- rbind(pop_data, df)
    }
  )
  
  # 3. Crear grilla para visualizar la función Rastrigin
  x <- seq(-5.12, 5.12, length.out = 100)
  y <- seq(-5.12, 5.12, length.out = 100)
  grid <- expand.grid(X1 = x, X2 = y)
  grid$Z <- apply(grid, 1, func)
  
  # 4. Graficar y animar
  base_plot <- ggplot() +
    geom_raster(data = grid, aes(x = X1, y = X2, fill = Z), interpolate = TRUE) +
    scale_fill_viridis_c() +
    geom_point(data = pop_data, aes(x = X1, y = X2), color = "red", size = 1, alpha = 0.6) +
    labs(title = "Optimización de Rastrigin usando GA", subtitle = "Generación: {closest_state}",
         x = "x1", y = "x2") +
    transition_states(Generacion, transition_length = 2, state_length = 1) +
    theme_minimal()
  
  # 5. Exportar como GIF
  anim_save("optim_rastrigin_ga.gif", animation = animate(base_plot, renderer = gifski_renderer(), fps = 5, width = 600, height = 500))
}

Analizaremos cada función en 2 y 3 dimensiones, graficando su paisaje antes de la optimización y luego aplicando un AG con múltiples corridas para evaluar la robustez de los resultados.

1.6.1 Implementación en R de algoritmos evolutivos

Se utiliza la función ga() del paquete GA. Para problemas de minimización se define la función de fitness como el negativo del valor objetivo, ya que ga() maximiza por defecto. Se especifican los límites de búsqueda.

Los resúmenes en cada ejecución reportan el mejor fitness encontrado (negativo) y la solución óptima en cada ejecución.

1.6.2 Optimización de la función de Rosenbrock en 2 dimensiones

# Ejecución del método
ga_ros2d <- ga(type = "real-valued",
               fitness = function(x) -f_rosenbrock(x),
               lower = c(-5, -5), upper = c(5, 5),
               popSize = 50, maxiter = 100, run = 50)
summary(ga_ros2d)

# Graficación del proceso de optimización
gif_name <- "optim_rosenbrock_ga.gif"
animate_ga_optimization(f_rosenbrock)

set.seed(123)  # semilla reproducible
best_vals_ros <- replicate(30, {
  GA <- ga(type = "real-valued",
           fitness = function(x) -f_rosenbrock(x),
           lower = c(-5, -5), upper = c(5, 5),
           popSize = 50, maxiter = 100, run = 50)
  -GA@fitnessValue  # convertir a valor positivo
})
mean_ros <- mean(best_vals_ros)
sd_ros   <- sd(best_vals_ros)

1.6.3 Optimización de la función de Rosenbrock en 3 dimensiones

# Ejecución del método
ga_ros3d <- ga(type = "real-valued",
               fitness = function(x) -f_rosenbrock(x),
               lower = c(-5, -5, 3), upper = c(5, 5, 3),
               popSize = 50, maxiter = 100, run = 50)
summary(ga_ros3d)

# Realizar 30 ejecuciones independientes para Rosenbrock 3D
set.seed(123)  # semilla reproducible
best_vals_ros3d <- replicate(30, {
  GA <- ga(type = "real-valued",
           fitness = function(x) -f_rosenbrock(x),
           lower = c(-5, -5, 3), upper = c(5, 5, 3),
           popSize = 50, maxiter = 100, run = 50)
  -GA@fitnessValue  # convertir a valor positivo
})
mean_ros3d <- mean(best_vals_ros3d)
sd_ros3d   <- sd(best_vals_ros3d)
#  Rosenbrock 3D, Rastrigin 3D.

1.6.4 Optimización de la función de Rastrigin en 2 dimensiones

# Ejecución del método
ga_ras2d <- ga(type = "real-valued",
               fitness = function(x) -f_rastrigin(x),
               lower = c(-5, -12), upper = c(5, 12),
               popSize = 50, maxiter = 100, run = 50)
summary(ga_ras2d)

# Graficación del proceso de optimización
gif_name <- "optim_rastrigin_ga.gif"
animate_ga_optimization(f_rastrigin)

set.seed(123)  # semilla reproducible
best_vals_ras2d <- replicate(30, {
  GA <- ga(type = "real-valued",
           fitness = function(x) -f_rastrigin(x),
           lower = c(-5, -5), upper = c(5, 5),
           popSize = 50, maxiter = 100, run = 50)
  -GA@fitnessValue  # convertir a valor positivo
})
mean_ras2d <- mean(best_vals_ras2d)
sd_ras2d   <- sd(best_vals_ras2d)

1.6.5 Optimización de la función de Rastrigin en 3 dimensiones

# Ejecución del método
ga_ras3d <- ga(type = "real-valued",
               fitness = function(x) -f_rastrigin(x),
               lower = c(-5, -12,3), upper = c(5, 12,3 ),
               popSize = 50, maxiter = 100, run = 50)
summary(ga_ras3d)

set.seed(123)  # semilla reproducible
best_vals_ras3d <- replicate(30, {
  GA <- ga(type = "real-valued",
           fitness = function(x) -f_rastrigin(x),
           lower = c(-5, -5, 3), upper = c(5, 5, 3),
           popSize = 50, maxiter = 100, run = 50)
  -GA@fitnessValue  # convertir a valor positivo
})
mean_ras3d <- mean(best_vals_ras3d)
sd_ras3d   <- sd(best_vals_ras3d)

1.6.6 Cálculo de estadísticas y análisis

Para evaluar la variabilidad del método estocástico, se repite cada caso al menos 30 veces con semillas distintas. Se registra el mejor valor de fitness (valorizado positivamente) obtenido en cada corrida. Se calcularán estadísticas (media y desviación estándar) del mejor valor de fitness obtenido en 30 ejecuciones independientes de cada caso y se resumirán en una tabla.

Con los vectores de mejores valores (best_vals_ros, etc.), se calculan la media y desviación estándar de cada conjunto de 30 resultados. Por ejemplo, mean_ros y sd_ros arriba y las demas, Para asi presentar los resultados.

library(knitr)
resultados <- data.frame(
  Función   = c("Rosenbrock", "Rastrigin", "Rosenbrock", "Rastrigin"),
  Dimensión = c("2D", "2D", "3D", "3D"),
  Media     = c(mean_ros, mean_ras2d, mean_ros3d, mean_ras3d),
  SD        = c(sd_ros, sd_ras2d, sd_ros3d, sd_ras3d)
)
kable(resultados, caption = "Resumen estadístico (media y desviación estándar) del mejor fitness obtenido tras 30 ejecuciones independientes de cada caso.")

Los resultados de las múltiples ejecuciones se resumen en la Tabla 1. Esta tabla muestra la media y desviación estándar del mejor valor de fitness (recordado que es el valor de la función objetivo en su mínimo global, típicamente cercano a 0) para cada combinación de función y dimensión. Se observa que para Rosenbrock 2D, la media del fitness mínimo es cercana a 0 con baja dispersión, reflejando que el GA normalmente encuentra el mínimo global (0) o cercano. Para Rastrigin 2D, la media también puede acercarse a 0, pero con mayor desviación estándar debido a los múltiples mínimos locales. En 3D ambos problemas suelen mostrar valores medios mayores (más alejados de 0) y mayor variabilidad, lo cual indica una mayor dificultad de búsqueda al aumentar la dimensionalidad.

# tabla de los valores calculados)
library(knitr)
res_df <- data.frame(
  Función   = c("Rosenbrock", "Rastrigin", "Rosenbrock", "Rastrigin"),
  Dimensión = c("2D", "2D", "3D", "3D"),
  Media     = c(mean_ros, mean_ras2d, mean_ros3d, mean_ras3d),
  SD        = c(sd_ros, sd_ras2d, sd_ros3d, sd_ras3d)
)
kable(res_df, caption = "Tabla 1. Estadísticas (media y desviación estándar) del fitness mínimo alcanzado en 30 corridas independientes para cada función y dimensión.")

Tabla 1. Estadísticas (media y desviación estándar) del fitness mínimo alcanzado en 30 corridas independientes para cada función y dimensión.

1.6.7 Conclusiones método de algoritmos evolutivos

Los resultados confirman que el algoritmo genético es capaz de aproximarse a los mínimos globales de ambos problemas en múltiples dimensiones. Como era de esperar, Rastrigin mostró mayor variabilidad en los valores de fitness debido a sus muchos mínimos locales, lo que implica que algunas ejecuciones del GA pueden quedarse atrapadas en óptimos locales alejados del global. En contraste, Rosenbrock (aunque es no convexa) tiende a un único valle principal; por ello, la mayoría de las corridas alcanzaron valores cercanos al mínimo global con menor dispersión. En general se observa que al aumentar la dimensión (de 2D a 3D) la tarea se complica y la media del fitness aumenta (peor óptimo encontrado), reflejando la maldición de la dimensionalidad. El uso de múltiples ejecuciones independientes es esencial para evaluar la robustez de los AG. Debido a su naturaleza estocástica, cada ejecución puede converger a soluciones distintas. Al analizar la media y desviación estándar de los fitness finales se obtiene una medida de fiabilidad del algoritmo: una baja desviación indica resultados consistentes. En la literatura sobre algoritmos genéticos se reconoce que en muchos casos una sola ejecución puede no ser representativa​jstatsoft.org. Aunque un análisis comparativo profundo (p.ej., usando poblaciones más grandes o múltiples corridas en paralelo) queda fuera del alcance de este documento, nuestros resultados ilustran este fenómeno. Este estudio es reproducible: todo el código R necesario está incluido, permitiendo a otros investigadores replicar los experimentos, variar parámetros del GA (tasa de cruce, mutación, tamaño de población, etc.) y comparar con otros algoritmos de optimización.:Conclusiones Se ha presentado una documentación completa de la optimización de las funciones de Rosenbrock y Rastrigin en 2D y 3D empleando algoritmos genéticos en R. Mediante visualizaciones 3D iniciales se ilustraron las características de cada función de prueba. Se implementó el paquete GA para resolver cada caso y se realizaron 30 ejecuciones independientes para evaluar la robustez. Los resultados muestran que el GA puede encontrar aproximaciones al mínimo global en ambos problemas, aunque la función Rastrigin (múltiples mínimos locales) presenta más variabilidad y dificultad, especialmente en 3D.

Parte 2. Optimización combinatoria

2.1 Planteamiento del problema

El problema se puede plantear como una instancia particular del problema del viajero en ciencias de la computación. En la Tabla 2 se muestran las 13 ciudades principales de Colombia con su respectivas latitudes y longitudes.

	Latitud	Longitud
Bogotá	4.7110	-74.0721
Medellín	6.2442	-75.5812
Cali	3.4516	-76.5320
Barranquilla	10.9685	-74.7813
Cartagena	10.3910	-75.4794
Cúcuta	7.8941	-72.5078
Soledad	10.9264	-74.8055
Ibagué	4.4389	-75.2322
Bucaramanga	7.1193	-73.1227
Villavicencio	4.1420	-73.6298
Santa Marta	11.2408	-74.1990
Manizales	5.0703	-75.5138
Pereira	4.8143	-75.6946

Tabla 2. Ciudades principales de Colombia junto con su latencia y longitud.

El objetivo de este problema es optimizar con respecto a los costos y no a las distancias entre ciudades, por lo que se tienen que considerar los costos de combustible, de peajes y cuánto cobra el vendedor por hora de trabajo. Sin embargo, la distancia afectará al valor de todos estos costos, y por lo tanto hay que considerarla. Adicionalmente, todos estos costos están en función de la distancia en metros, por lo que es necesario realizar la conversión de latitud-longitud a coordenadas en metros. Para esto se usó la librería “sf” para realizar la conversión al sistema de coordenadas planas UTM, el cual representa los puntos del globo en metros.

La tabla de coordenadas en metros utilizando la información de la Tabla 1 se implementa a continuación:

nombre_ciudades <- c("Bogotá", "Medellín", "Cali", "Barranquilla", "Cartagena", 
             "Cúcuta", "Pasto", "Ibagué", "Bucaramanga", "Villavicencio", 
             "Santa Marta", "Manizales", "Pereira")

# Definición de ciudades, latitud y longitud de cada ciudad.
ciudades <- tibble::tibble(
  Ciudad = nombre_ciudades,
  Latitud = c(4.7110, 6.2442, 3.4516, 10.9685, 10.3910, 
              7.8941, 1.2136, 4.4389, 7.1193, 4.1420, 
              11.2408, 5.0703, 4.8143),
  Longitud = c(-74.0721, -75.5812, -76.5320, -74.7813, -75.4794, 
               -72.5078, -77.2811, -75.2322, -73.1227, -73.6298, 
               -74.1990, -75.5138, -75.6946)
)

# Time zone value
# Por simplicidad, se elige el mismo punto común para calcular el timezone.
# Se utiliza la longitud de Bogotá como punto común para calcular el utm_zone.
utm_zone <- floor((-74.0721 + 180)/6) + 1 + 32700

# Convertimos a objeto espacial
puntos <- st_as_sf(ciudades, coords = c("Longitud", "Latitud"), crs = 4326)

puntos_en_utm <- st_transform(puntos, crs=utm_zone)
ciudades_en_metros <- st_coordinates(puntos_en_utm)
tabla_coordenadas_ciudades <- data.frame(
  ciudad = nombre_ciudades,
  coord_x_en_metros = ciudades_en_metros[,1],
  coord_y_en_metros = ciudades_en_metros[,2]
)
tabla_coordenadas_ciudades

Luego, se crea la matriz de distancias en metros entre cada ciudad.

coordenadas_ciudades <- data.matrix(tabla_coordenadas_ciudades[,c(2,3)])
distancias <- compute_distance_matrix(coordenadas_ciudades)
print(distancias)

Ahora, hay que considerar los costos para poder crear una matriz de costos para resolver el problema de optimización.

Con respecto al costo del combustible, este dependerá del vehículo que el vendedor utilizará para realizar las entregas. Para este problema, se utilizará un furgón DFSK C35, el cual tiene un consumo de combustible de 7.6 litros aproximado por cada 100 Km de recorrido, que se puede traducir en 0.000076 litros por cada 1 metro (DFSK, s.f., Especificaciones Furgon C35). Adicionalmente, se considera el precio promedio de la gasolina en Colombia que, segun (La República, 2025), tiene un valor de $15.827 pesos colombianos por galón, que se puede traducir en $4.022 pesos colombianos por litro aproximadamente.

gasto_litro_por_metro <- 0.000076
precio_gasolina_por_litro <- 4.022 
costos_combustible <- distancias * gasto_litro_por_metro * precio_gasolina_por_litro
rownames(costos_combustible) <- nombre_ciudades
colnames(costos_combustible) <- nombre_ciudades
print(costos_combustible)

Considerando el costo de los peajes, se usarán los datos de (Autofact, 2024) para calcular el valor promedio de cada peaje de las ciudades consideradas y se asumirá que el vendedor pagará el valor promedio del peaje de la ciudad de origen y de la ciudad de destino exactamente 1 vez cada que las recorra.

matriz_costos <- costos_vendedor + costos_peajes + costos_combustible
print(matriz_costos)

                 Bogotá  Medellín      Cali Barranquilla Cartagena    Cúcuta     Pasto    Ibagué Bucaramanga Villavicencio Santa Marta Manizales   Pereira
Bogotá             0.00  67240.80  74529.48    137649.07 131493.79  91153.12 113638.21  48319.61    73468.06      37405.20   143474.36  54767.84  58441.36
Medellín       67240.80      0.00  80586.62    112936.95 103200.50  92810.96 126574.14  62860.46    76622.06      79451.17   121635.15  51827.65  57648.70
Cali           74529.48  80586.62      0.00    162359.33 151607.74 134676.25  68676.65  55156.30   116491.41      77704.02   171434.62  61299.49  56681.47
Barranquilla  137649.07 112936.95 162359.33         0.00  38637.58  93566.43 208193.11 143529.48   100192.54     148166.46    33370.49 133613.92 139789.55
Cartagena     131493.79 103200.50 151607.74     38637.58      0.00  96527.86 197108.70 134796.68    99220.12     142737.87    51591.87 124251.54 130100.20
Cúcuta         91153.12  92810.96 134676.25     93566.43  96527.86      0.00 178084.43 108103.99    45478.74      97007.68    93844.11 104121.58 110845.74
Pasto         113638.21 126574.14  68676.65    208193.11 197108.70 178084.43      0.00  98011.17   159865.59     111610.69   217401.50 106738.98 102216.29
Ibagué         48319.61  62860.46  55156.30    143529.48 134796.68 108103.99  98011.17      0.00    89865.26      55455.92   151187.62  41556.84  40924.53
Bucaramanga    73468.06  76622.06 116491.41    100192.54  99220.12  45478.74 159865.59  89865.26        0.00      80445.13   103039.49  86147.79  92820.63
Villavicencio  37405.20  79451.17  77704.02    148166.46 142737.87  97007.68 111610.69  55455.92    80445.13          0.00   153360.70  65113.75  67537.40
Santa Marta   143474.36 121635.15 171434.62     33370.49  51591.87  93844.11 217401.50 151187.62   103039.49     153360.70        0.00 141825.81 148203.21
Manizales      54767.84  51827.65  61299.49    133613.92 124251.54 104121.58 106738.98  41556.84    86147.79      65113.75   141825.81      0.00  36827.52
Pereira        58441.36  57648.70  56681.47    139789.55 130100.20 110845.74 102216.29  40924.53    92820.63      67537.40   148203.21  36827.52      0.00

Por último, se debe de considerar el salario por horas promedio de un vendedor que haga dichos recorridos para entregar los pedidos. Adicional a ser vendedor, esta persona realizará las entregas y la conducción de la mercancía, por lo que también hay que considerar su pago como transportista. Dicho esto, se promediaron salarios de un transportista y un vendedor para poder determinar el salario final de la persona. Según (Talent.com, s.f., Salario medio para Transporte De Carga en Colombia 2025) y (Talent.com, s.f., Salario medio para Vendedor en Colombia 2025), un transportista gana en promedio $9.341 pesos colombianos la hora y un vendedor hace $7.144 la hora; por lo tanto, el salario a utilizar será el $8.242 pesos la hora.

Asumiendo que el vendedor conducirá a una velocidad media de 50 km/h al día, que se podrían traducir en 50000 m/h para facilitar cálculos, se procede a calcular el tiempo de viaje y, junto con esto, cuánto se le pagará al vendedor por cada viaje.

velocidad_media_metros_por_hora <- 50000
salario_vendedor_colombia <- 7144
salario_transportista_colombia <- 9341
salario_final <- (salario_vendedor_colombia + salario_transportista_colombia)/2
costos_vendedor <- distancias*(1/velocidad_media_metros_por_hora)*salario_final
rownames(costos_vendedor) <- nombre_ciudades
colnames(costos_vendedor) <- nombre_ciudades
costos_vendedor

Concluyendo, la matriz de costos simplemente sería el resultado de la suma de las anteriores matrices previamente definidas, y esta será la matriz que se utilizarán en los distintos algoritmos para encontrar la mejor ruta.

matriz_costos <- costos_vendedor + costos_peajes + costos_combustible
print(matriz_costos)

2.2 Implementación Método de Colonia de hormigas

A continuación, se implementa el método de colonia de hormigas para encontrar

# Número de ciudades
n_ciudades <- 13

# Función para normalizar la matriz de costos, ya que search_tour_ants
# funciona con valores normalizados.
# Se optó por usar normalización zscore.
normalize_zscore <- function(m) {
  (m - mean(m)) / sd(m)
}

# Normalización de matriz de costos
matriz_costos_normalizada <- normalize_zscore(matriz_costos)

# Ejecución del método de colonia de hormigas
recorrido_optimizado <- search_tour_ants(matriz_costos_normalizada, n_ciudades, K = 80, N = 40, log=TRUE)
print("MEJOR RECORRIDO: ")
print(recorrido_optimizado$tour)

# Tour más óptimo: 1 10  7  3  8 13 12  2  5  4 11  6  9

# Rutas seleccionadas entre algunas ciudades
print("RUTA IDEAL: ")
print("Bogotá->Villavicencio->Pasto->Cali->Ibagué->Pereira->Manizales->Medellín->Cartagena->Barranquilla->Santa Marta->Cúcuta->Bucaramanga")
pares_rutas <- list(
  c("Bogotá", "Villavicencio"),
  c("Villavicencio", "Pasto"),
  c("Pasto", "Cali"),
  c("Cali", "Ibagué"),
  c("Ibagué", "Pereira"),
  c("Pereira", "Manizales"),
  c("Manizales", "Medellín"),
  c("Medellín", "Cartagena"),
  c("Cartagena", "Barranquilla"),
  c("Barranquilla", "Santa Marta"),
  c("Santa Marta", "Cúcuta"),
  c("Cúcuta", "Bucaramanga")
)

# Convertimos a objeto espacial
puntos <- st_as_sf(ciudades, coords = c("Longitud", "Latitud"), crs = 4326)

# Obtener rutas con osrm
rutas <- list()
for (par in pares_rutas) {
  origen <- puntos[ciudades$Ciudad == par[1], ]
  destino <- puntos[ciudades$Ciudad == par[2], ]
  ruta <- try(osrmRoute(src = origen, dst = destino, ), silent = TRUE)
  if (!inherits(ruta, "try-error")) {
    rutas <- append(rutas, list(ruta))
  }
}

# Crear el mapa
mapa <- leaflet() %>%
  addTiles() %>%
  addCircleMarkers(data = puntos, label = ~Ciudad, radius = 6, color = "blue", fillOpacity = 0.8)

# Añadir las rutas en rojo
for (ruta in rutas) {
  mapa <- mapa %>% addPolylines(data = ruta, color = "red", weight = 3)
}

print("COSTO FINAL DEL RECORRIDO: ")
costos <- c(
  matriz_costos[1,10], #"Bogotá" -> "Villavicencio"
  matriz_costos[10,7], #"Villavicencio" -> "Pasto"
  matriz_costos[7,3], #"Pasto" -> "Cali"
  matriz_costos[3,8], #"Cali" -> "Ibagué"
  matriz_costos[8,13], #"Ibagué" -> "Pereira"
  matriz_costos[13,12], #"Pereira" -> "Manizales"
  matriz_costos[12,2], #"Manizales" -> "Medellín"
  matriz_costos[2,5], #"Medellín" -> "Cartagena"
  matriz_costos[5,4], #"Cartagena" -> "Barranquilla"
  matriz_costos[4,11], #"Barranquilla" -> "Santa Marta"
  matriz_costos[11,6], #"Santa Marta" -> "Cúcuta"
  matriz_costos[6,9] #"Cúcuta" -> "Bucaramanga"
)
print(sum(costos))
mapa

2.3 Conclusiones

Conclusiones finales

Reporte de contribución individual

Leonardo Federico Corona Torres

• Apoyo en redacción y realización de implementaciones de sección “Método de algoritmos evolutivos”.

• Apoyo en redacción de implementaciones de graficación en Parte 2 del reporte.

• Búsqueda e investigación de funciones de prueba a utilizar.

• Apoyo en redacción y realización de implementaciones de sección “Selección e implementación de las funciones de prueba”.

David Escobar Ruiz

• Definición y escritura de estructura del reporte.

• Apoyo en redacción y realización de implementaciones de sección “Selección e implementación de las funciones de prueba”.

• Apoyo en redaccióny realización de implementaciones de sección “Método de optimización de partículas”.

• Apoyo en redacción y realización de implementaciones de métodos, planteamiento y preprocesamiento de datos Parte 2 del reporte.

Johan Sebastián Robles Rincón

• Implementación y prueba del método de evolución diferencial en R para funciones de Rosenbrock y Rastrigin en 2D y 3D.

• Análisis de resultados numéricos y gráficos.

• Apoyo en redacción de sección del contenido técnico relacionado con este método de evolucion diferencial en el informe.

Sebastian Soto Arcila

• Apoyo en redacción y realización de implementaciones de sección “Método de optimización por descenso de gradiente”

• Apoyo en redacción de sección “Selección e implementación de las funciones de prueba”.

Repositorio de GitHub del proyecto

https://github.com/druiz35/RNABI2025-1-Equipo3/

Bibliografía

Autofact. (2024). Peajes en Colombia: Conoce sus ubicaciones y precios 2024. https://www.autofact.com.co/blog/mi-carro/peajes/peajes-colombia-precios

DFSK. (s.f.). Especificaciones Furgon C35. Recuperado el 2 de mayo de 2025, de http://static.multiaviso.com/vehicle/specs/22-VRZC564XLBE6-dfsk-otros-modelos-2017-furgon-c35.pdf

La República. (2025). PRECIO DE LA GASOLINA. Recuperado el 2 de mayo de 2025, de https://www.larepublica.co/precio-de-la-gasolina

Molga, Smutnicki. (2005). Test functions for optimization needs. https://robertmarks.org/Classes/ENGR5358/Papers/functions.pdf

Talent.com. (s.f.). Salario medio para Transporte De Carga en Colombia 2025. Recuperado el 2 de mayo de 2025, de https://co.talent.com/salary?job=transporte+de+carga

Talent.com. (s.f.). Salario medio para Vendedor en Colombia 2025. Recuperado el 2 de mayo de 2025, de https://co.talent.com/salary?job=vendedor

Wikipedia. (s.f.). Test functions for optimization. Wikipedia. Recuperado el 2 de mayo de 2025, de https://en.wikipedia.org/wiki/Test_functions_for_optimization

Wikipedia. (s.f.). Rosenbrock function. Wikipedia. Recuperado el 2 de mayo de 2025, de https://en.wikipedia.org/wiki/Rosenbrock_

Wikipedia. (s.f.). Rastrigin function. Wikipedia. Recuperado el 2 de mayo de 2025, de https://en.wikipedia.org/wiki/Rastrigin_function

X.-S. Yang, Test problems in optimization, in: Engineering Optimization: An Introduction with Metaheuristic Applications (Eds Xin-She Yang), John Wiley & Sons, (2010)