Parte 1: Optimización Numérica

Optimice las funciones en dos y tres dimensiones usando un método de descenso por gradiente con condición inicial aleatoria
Optimice las funciones en dos y tres dimensiones usando: algoritmos evolutivos, optimización de partículas y evolución diferencial
Represente con un gif animado o un video el proceso de optimización de descenso por gradiente y el proceso usando el método heurístico.
¿Qué aportaron los métodos de descenso por gradiente y qué aportaron los métodos heurísticos? Para responder a esta pregunta considere el valor final de la función objetivo y el número de evaluaciones de la función objetivo. Para responder a esta pregunta es posible que se requiera hacer varias corridas de los algoritmos.

0. Exploración de las funciones

Función de Schwefel

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

Descripción:

Dimensiones: \(d\)

La función de Schwefel cuenta con múltiples mínimos locales.

Dominio:

La función se evalúa en el hipercubo: \(x_i \in [-500, 500]\), para todo \(i = 1, \dots, d\).

Mínimo Global:

\(f(\mathbf{x}^*) = 0\), en \(\mathbf{x}^* = (420.9687, \dots, 420.9687)\)

Definición de Función

# Función Schwefel
schwefel <- function(x) {
  A <- 418.9829
  d <- length(x)
  z <- A*d - sum(x * sin(sqrt(abs(x))))
 
  return (z)
}

Sabemos que la función de Schwefel tiene un mínimo local en \(x^*=(420.9687, \cdots, 420.9687)\)

schwefel(c(420.9687, 420.9687))

## [1] 2.545567e-05

schwefel(c(420.9687, 420.9687, 420.9687))

## [1] 3.818351e-05

Estos resultados representan valores aproximados a cero debido a errores de redondeo introducidos en la programación de la función.

Ahora vamos a visualizar la función:

# Código para graficar cualquier función
plot_function <- function(func, x_range = c(-500, 500), n_points = 1000, highlight_x = NULL) {
  # Crear el vector de x en el rango especificado
  x <- seq(x_range[1], x_range[2], length.out = n_points)
 
  # Evaluar la función en los puntos de x
  y <- sapply(x, function(xi) func(c(xi)))
 
  # Graficar la función
  plot(x, y, type = "l", col = "blue", lwd = 2,
       main = paste("Gráfico de la función", deparse(substitute(func))),
       xlab = "x", ylab = paste(deparse(substitute(func)), "(x)"))
 
  # Si se especifica un valor de x para destacar, agregar la línea vertical y el punto
  if (!is.null(highlight_x)) {
    highlight_y <- func(c(highlight_x))
    abline(v = highlight_x, col = "red", lty = 2)
    points(highlight_x, highlight_y, col = "red", pch = 19)
  }
}

plot_function(schwefel, highlight_x = 420.9687)

Fig 1. Evaluación de Schwefel (2D)

plot_function_3d <- function(func,
                           x_range = c(-500, 500),
                           y_range = c(-500, 500),
                           resolution = 100,
                           optimal_points = NULL,
                           colorscale = "Viridis",
                           title = "Visualización 3D de función",
                           x_label = "x",
                           y_label = "y",
                           z_label = "f(x,y)") {
 
  # Crear grid de puntos
  x <- seq(x_range[1], x_range[2], length.out = resolution)
  y <- seq(y_range[1], y_range[2], length.out = resolution)
 
  # Evaluar la función en toda la cuadrícula
  z <- outer(x, y, Vectorize(function(x, y) func(c(x, y))))
 
  # Crear el gráfico 3D interactivo
  p <- plot_ly(x = ~x, y = ~y, z = ~z) %>%
    add_surface(colorscale = colorscale)
 
  # Agregar puntos óptimos si se especifican
  if (!is.null(optimal_points)) {
    for (i in 1:nrow(optimal_points)) {
      point <- optimal_points[i, ]
      p <- p %>% add_markers(
        x = point$x,
        y = point$y,
        z = func(c(point$x, point$y)),
        marker = list(color = point$color, size = point$size),
        name = point$name
      )
    }
  }
 
  # Configurar el diseño del gráfico
  p <- p %>% layout(
    title = title,
    scene = list(
      xaxis = list(title = x_label),
      yaxis = list(title = y_label),
      zaxis = list(title = z_label)
    )
  )
 
  return(p)
}

# Definir óptimo schwefel
optimal_points <- data.frame(
  x = c(420.9687),
  y = c(420.9687),
  name = c("Óptimo global"),
  color = c("red"),
  size = c(5)
)

# Graficar la función con su óptimo
plot_function_3d(
  func = schwefel,
  x_range = c(-500, 500),
  y_range = c(-500, 500),
  resolution = 100,
  optimal_points = optimal_points,
  title = "Función de Schwefel (2D -> 3D)",
  z_label = "Schwefel(x, y)"
)

Fig 2. Evaluación de Schwefel (3D)

# Crear una grilla de puntos (para 2D)
x_vals <- seq(-500, 500, length.out = 200)
y_vals <- seq(-500, 500, length.out = 200)
grid <- expand.grid(x = x_vals, y = y_vals) %>%
  rowwise() %>%
  mutate(z = schwefel(c(x, y))) %>%
  ungroup()

# Crear gráfico de curvas de nivel
ggplot(grid, aes(x = x, y = y, z = z)) +
  geom_contour_filled(bins = 30) +  # Colores suaves
  geom_contour(color = "black", alpha = 0.3) +  # Líneas de nivel
  geom_point(aes(x = 420.9687, y = 420.9687), color = "red", size = 3) +  # Óptimo
  annotate("text", x = 420.9687, y = 420.9687 + 30, label = "Óptimo", color = "red") +
  theme_minimal() +
  scale_fill_viridis_d() +
  labs(
    title = "Curvas de nivel - Función de Schwefel",
    subtitle = "Visualización 2D de la superficie 3D",
    x = "Valor de x_1",
    y = "Valor de x_2",
    fill = "Valor de la función"
  )

Fig 3. Curvas de Nivel de Schwefel

Función de Griewank

\[ f(\mathbf{x}) = 1 + \frac{1}{4000} \sum_{i=1}^{d} x_i^2 - \prod_{i=1}^{d} \cos\left(\frac{x_i}{\sqrt{i}}\right) \]

Descripción:

Dimensiones: \(d\)

La función de Griewank tiene una estructura ondulada con múltiples mínimos locales, pero un único mínimo global en el origen.

Dominio:

\(x_i \in [-600, 600]\), para todo \(i = 1, \dots, d\)

Mínimo Global:

\(f(\mathbf{0}) = 0\), en \(\mathbf{x} = (0, \dots, 0)\)

Definición de Función

# Función Griewank (2D y 3D)
griewank <- function(x) {
  sum_term <- sum(x^2) / 4000
  prod_term <- prod(cos(x / sqrt(seq_along(x))))
  return(1 + sum_term - prod_term)
}

griewank(c(0, 0))

## [1] 0

griewank(c(0, 0, 0))

## [1] 0

Estos resultados representan el mínimo global de la función de Griewank, que es 0 en el origen.

plot_function(griewank, x_range = c(-600, 600), highlight_x = 0)

# Definir óptimo griewank
optimal_points <- data.frame(
  x = c(0),
  y = c(0),
  name = c("Óptimo global"),
  color = c("red"),
  size = c(5)
)

# Graficar la función con su óptimo
plot_function_3d(
  func = griewank,
  x_range = c(-600, 600),
  y_range = c(-600, 600),
  resolution = 100,
  optimal_points = optimal_points,
  title = "Función de Griewank (2D -> 3D)",
  z_label = "Griewank(x, y)"
)

Fig 4. Evaluación de Griewank (3D)

# Crear una grilla de puntos (para 2D)
x_vals <- seq(-600, 600, length.out = 200)
y_vals <- seq(-600, 600, length.out = 200)
grid <- expand.grid(x = x_vals, y = y_vals) %>%
  rowwise() %>%
  mutate(z = griewank(c(x, y))) %>%
  ungroup()

# Crear gráfico de curvas de nivel
ggplot(grid, aes(x = x, y = y, z = z)) +
  geom_contour_filled(bins = 30) +  # Colores suaves
  geom_contour(color = "black", alpha = 0.3) +  # Líneas de nivel
  geom_point(aes(x = 0, y = 0), color = "red", size = 3) +  # Óptimo
  annotate("text", x = 0, y = 30, label = "Óptimo", color = "red") +
  theme_minimal() +
  scale_fill_viridis_d() +
  labs(
    title = "Curvas de nivel - Función de Griewank",
    subtitle = "Visualización 2D de la superficie 3D",
    x = "Valor de x_1",
    y = "Valor de x_2",
    fill = "Valor de la función"
  )

Fig 5. Curvas de Nivel de Griewank

1. Optimización por descenso del gradiente o métodos Quasi-Newton

Como podemos ver, la función cuenta con varios óptimos locales que pueden hacer que un método numérico basado en el gradiente sea propenso a errores. Para ello vamos a ver cómo se comporta la solución del método optim(method = "BFGS") que es el más parecido al que vimos en clase. Este método Quasi-Newton optimiza la función utilizando una aproximación del gradiente y la matriz Hessiana, sin calcularlos directamente.

Antes de esto vamos a crear vectores para seguir las soluciones halladas.

# Vectores para almacenar los valores intermedios de la función
schwefel_2d_history <- list()

# Función para llenar el camino de valores 2D
schwefel_with_history_2d <- function(x) {
  val <- schwefel(x)
  schwefel_2d_history <<- c(schwefel_2d_history, list(list(x = x, value = val)))  # Guarda el valor intermedio
  return(val)
}

schwefel_3d_history <- list()

# Función para llenar el camino de valores 3D
schwefel_with_history_3d <- function(x) {
  val <- schwefel(x)
  schwefel_3d_history <<- c(schwefel_3d_history, list(list(x = x, value = val)))  # Guarda el valor intermedio
  return(val)
}

# Optimización 2D

set.seed(1987)

x0_2d <- runif(1, min=-500, max=500)

# Optimización
res_2d <- optim(
  par = x0_2d,
  fn = schwefel_with_history_2d,
  method = "BFGS",
  control = list(trace = 1, REPORT = 1, maxit = 500),
  hessian = FALSE  # En 2d no usamos el Hessiano sino el gradiente únicamente
)

## initial  value 790.168681 
## iter   2 value 765.448307
## iter   3 value 729.639629
## iter   4 value 679.713597
## iter   5 value 613.806275
## iter   6 value 533.106774
## iter   7 value 443.514431
## iter   8 value 355.131538
## iter   9 value 138.514178
## iter  10 value 122.084251
## iter  11 value 120.486081
## iter  12 value 118.438625
## iter  13 value 118.438347
## iter  13 value 118.438347
## final  value 118.438347 
## converged

# Resultado
res_2d$par       # Coordenadas del mínimo encontrado

## [1] -302.5247

res_2d$value     # Valor mínimo de la función

## [1] 118.4383

# Optimización 3D

set.seed(1987)

x0_3d <- runif(2, min=-500, max=500)

# Optimización
res_3d <- optim(
  par = x0_3d,
  fn = schwefel_with_history_3d,
  method = "BFGS",
  control = list(trace = 1, REPORT = 1, maxit = 1000),
)

## initial  value 1508.888943 
## iter   2 value 1483.707203
## iter   3 value 1447.173577
## iter   4 value 1396.108891
## iter   5 value 1328.414538
## iter   6 value 1244.914455
## iter   7 value 1150.944012
## iter   8 value 1055.745413
## iter   9 value 473.273255
## iter  10 value 398.813939
## iter  11 value 398.132357
## iter  12 value 394.929298
## iter  13 value 394.916344
## iter  14 value 394.899954
## iter  14 value 394.899953
## iter  14 value 394.899953
## final  value 394.899953 
## converged

# Resultado
res_3d$par       # Coordenadas del mínimo encontrado

## [1] -25.8788 420.9689

res_3d$value     # Valor mínimo de la función

## [1] 394.9

2. Métodos Heurísticos

Primero definamos una lista que va a servir como el historial de cada método que nos va a servir posteriormente para hacer el gif.

2.1 Algoritmo Genético (GA)

# Optimización GA con historial (función genérica)
genetic_optimizer_with_history <- function(func, lower, upper, dim, popSize = 50, maxiter = 100) {
  history <- list()

  monitor <- function(obj) {
    best <- obj@population[which.max(obj@fitness), ]
    val <- func(best)
    history[[obj@iter]] <<- list(x = best, value = val)
  }

  result <- ga(
    type = "real-valued",
    fitness = function(x) -func(x),
    lower = rep(lower, dim),
    upper = rep(upper, dim),
    popSize = popSize,
    maxiter = maxiter,
    monitor = monitor
  )

  return(list(result = result, history = history))
}

ga_schwefel_2d <- genetic_optimizer_with_history(schwefel, -500, 500, dim = 2)
print(-ga_schwefel_2d$result@fitnessValue)

## [1] 8.197191

print(ga_schwefel_2d$result@solution)

##            x1       x2
## [1,] 418.6601 413.2229

ga_griewank_2d <- genetic_optimizer_with_history(griewank, -600, 600, dim = 2)
print(-ga_griewank_2d$result@fitnessValue)

## [1] 0.0003544587

print(ga_griewank_2d$result@solution)

##                 x1         x2
## [1,] -0.0008010589 -0.0376194

2.2 Optimización por Enjambre (PSO)

# PSO con historial
pso_optimizer_with_history <- function(func, lower, upper, dim, maxit = 100, s = 40) {
  history <- list()

  wrapped_func <- function(x) {
    val <- func(x)
    history[[length(history) + 1]] <<- list(x = x, value = val)
    return(val)
  }

  result <- psoptim(
    par = runif(dim, lower, upper),
    fn = wrapped_func,
    lower = rep(lower, dim),
    upper = rep(upper, dim),
    control = list(maxit = maxit, s = s)
  )

  return(list(result = result, history = history))
}

pso_schwefel_2d <- pso_optimizer_with_history(schwefel, -500, 500, dim = 2)
print(pso_schwefel_2d$result$value)

## [1] 2.545636e-05

print(pso_schwefel_2d$result$par)

## [1] 420.9687 420.9688

pso_griewank_2d <- pso_optimizer_with_history(griewank, -600, 600, dim = 2)
print(pso_griewank_2d$result$value)

## [1] 0.008588485

print(pso_griewank_2d$result$par)

## [1]  3.140275 -4.369408

2.3 Evolución Diferencial (DE)

3. Visualización de Resultados

3.1 Animación del Descenso de Gradiente

Para crear la animación debemos generar múltiples imágenes usando el historial de soluciones evaluadas en la función objetivo.

Schwefel

# Base de la curva de la función schwefel
x_vals <- seq(-500, 500, length.out = 1000)
curve_df_2d <- data.frame(x = x_vals, value = sapply(x_vals, function(xi) schwefel(c(xi))))

# Extraer historia del optimizador
schwefel_df_2d <- data.frame(
  step = 1:length(schwefel_2d_history),
  x = sapply(schwefel_2d_history, function(p) p$x),
  value = sapply(schwefel_2d_history, function(p) p$value)
)

Camino de la solución en Schwefel 2D

p <- ggplot() +
  geom_line(data = curve_df_2d, aes(x = x, y = value), color = "blue", linewidth = 1.2) +
  geom_point(data = schwefel_df_2d, aes(x = x, y = value), color = "red", size = 2) +
  transition_reveal(along = step) +
  labs(title = "Paso del optimizador: {frame_along}", x = "x", y = "Schwefel(x)") +
  theme_minimal()

animate(p, fps = 4, width = 600, height = 400, renderer = gifski_renderer("files/schwefel_path.gif"))

knitr::include_graphics("files/schwefel_path.gif")

Como podemos notar, los algoritmos numéricos dependen altamente de la elección del punto inicial \(x_0^*\). En el GIF se ve además que la solución se desplaza por fuera de la gráfica, esto se debe a que el optimizador podría estar explorando soluciones que se mueven en \(c(x1)\), pero al graficar los frames se pueden ver puntos intermedios en la solución. Sin embargo finalmente se vuelve a la gráfica cuando se está estabilizando el óptimo.

Camino de la solución en Schwefel 3D con curvas de nivel

x_vals <- seq(-500, 500, length.out = 200)
y_vals <- seq(-500, 500, length.out = 200)
grid <- expand.grid(x = x_vals, y = y_vals) %>%
  rowwise() %>%
  mutate(z = schwefel(c(x, y))) %>%
  ungroup()

schwefel_df_3d <- data.frame(
  step = 1:length(schwefel_3d_history),
  x = sapply(schwefel_3d_history, function(p) p$x[1]),
  y = sapply(schwefel_3d_history, function(p) p$x[2]),
  z = sapply(schwefel_3d_history, function(p) p$value)
)

p <- ggplot() +
  # Capa BASE de relleno
  geom_raster(data = grid, aes(x = x, y = y, fill = z), alpha = 0.5) +
  scale_fill_viridis_c(option = "viridis", name = "Valor Schwefel") +
  
  # Curvas de nivel
  geom_contour(data = grid, aes(x = x, y = y, z = z), 
               color = "black", alpha = 0.4, bins = 15) +
  
  # Trayectoria del optimizador
  geom_path(data = schwefel_df_3d, aes(x = x, y = y), 
            color = "darkgrey", alpha = 0.5, linewidth = 0.8) +
  
  # Puntos animados (nuevo esquema de colores)
  geom_point(data = schwefel_df_3d, aes(x = x, y = y, color = step), 
             size = 3, show.legend = FALSE) +
  scale_color_gradient(low = "orange", high = "red") +  # Azul a verde
  
  # Punto óptimo
  geom_point(aes(x = 420.9687, y = 420.9687), 
             shape = 17, color = "cyan", size = 4) +
  annotate("text", x = 420.9687, y = 450, label = "Óptimo Global", 
           color = "cyan", fontface = "bold") +
  
  # Estilo y animación
  theme_minimal() +
  labs(
    title = "Optimización de la Función Schwefel\nPaso: {frame_time}",
    x = "Coordenada x₁",
    y = "Coordenada x₂"
  ) +
  transition_time(step) +
  shadow_wake(wake_length = 0.1, alpha = 0.3) +
  theme(plot.margin = margin(t = 30, r = 10, b = 10, l = 10, unit = "pt"))

# Renderizar la animación
animate(p, fps = 7, width = 800, height = 600, renderer = gifski_renderer("files/schwefel_3d_path.gif"))

knitr::include_graphics("files/schwefel_3d_path.gif")

# COMPLETAR:
# 1. Crear dataframe con trayectorias
# 2. Graficar con transition_reveal()

3.2 Animación de Métodos Heurísticos

# Convertir historial en data.frame
convert_history_to_df <- function(history) {
  df <- data.frame(
    step = seq_along(history),
    x = sapply(history, function(p) p$x[1]),
    y = sapply(history, function(p) p$x[2]),
    value = sapply(history, function(p) p$value)
  )
  return(df)
}

# Base para curvas de nivel
create_grid <- function(func, lower, upper, n = 200) {
  x_vals <- seq(lower, upper, length.out = n)
  y_vals <- seq(lower, upper, length.out = n)
  grid <- expand.grid(x = x_vals, y = y_vals)
  grid$z <- apply(grid, 1, function(row) func(c(row[1], row[2])))
  return(grid)
}

# Crear y guardar animación
plot_trajectory_animation <- function(history_df, grid_df, title, output_file) {
  p <- ggplot() +
    geom_raster(data = grid_df, aes(x = x, y = y, fill = z), alpha = 0.5) +
    scale_fill_viridis_c() +
    geom_contour(data = grid_df, aes(x = x, y = y, z = z), color = "black", alpha = 0.4) +
    geom_path(data = history_df, aes(x = x, y = y), color = "darkgrey", alpha = 0.5) +
    geom_point(data = history_df, aes(x = x, y = y, color = step), size = 3, show.legend = FALSE) +
    scale_color_gradient(low = "orange", high = "red") +
    theme_minimal() +
    labs(title = title, x = "x", y = "y") +
    transition_time(step) +
    shadow_wake(wake_length = 0.1, alpha = 0.3)

  animate(p, fps = 6, width = 700, height = 500, renderer = gifski_renderer(output_file))
}

# Ejemplo de uso
history_df <- convert_history_to_df(ga_griewank_2d$history)
grid_df <- create_grid(griewank, -600, 600)
plot_trajectory_animation(history_df, grid_df, "GA Griewank 2D", "ga_griewank_2d.gif")

4. Conclusiones

# DISCUTIR:
# 1. Ventajas/desventajas de cada método
# 2. Sensibilidad a parámetros
# 3. Recomendaciones para cada función

TRABAJO 1: OPTIMIZACIÓN HEURÍSTICA

REDES NEURONALES Y ALGORITMOS BIOINSPIRADOS

Verificación de la carga de paquetes

Parte 1: Optimización Numérica

0. Exploración de las funciones

Función de Schwefel

Definición de Función

Función de Griewank

Definición de Función

1. Optimización por descenso del gradiente o métodos Quasi-Newton

2. Métodos Heurísticos

2.1 Algoritmo Genético (GA)

2.2 Optimización por Enjambre (PSO)

2.3 Evolución Diferencial (DE)

3. Visualización de Resultados

3.1 Animación del Descenso de Gradiente

Schwefel

3.2 Animación de Métodos Heurísticos

4. Conclusiones

Parte 2: Optimización Combinatoria

Reporte de Contribución Individual

Repositorio en Github

Bibliografía