Simulated Annealing

Introducción.

El Recocido Simulado (o Simulated Annealing) es un algoritmo de optimización inspirado en la metalurgia. Simula el comportamiento de un herrero tratando de conseguir la estructura perfecta en una espada: calienta el metal y luego lo enfría lentamente para que los átomos se acomoden de forma óptima.

En computación, se realiza lo mismo para encontrar la “mejor respuesta” a un problema complejo sin quedarse atrapados en malas soluciones.

¿Por qué calentar el problema?

Si se intenta resolver un problema yendo siempre hacia la mejor opción inmediata (esto se llama algoritmo constructivo, Greedy o “codicioso”), es muy probable que la solución se encuentre atrapado en un mínimo local. Es como intentar subir al Everest pero quedarse atrapado en la cima de una colina pequeña porque, para llegar a la montaña de verdad, primero se tendría que bajar un poco.

El Recocido Simulado resuelve esto permitiendo movimientos “malos” al principio.

Funcionamiento del proceso

El algoritmo sigue estos pasos fundamentales:

  1. Estado Inicial: se comienza con una solución aleatoria (o previa \(S\)) y una temperatura (\(T\)) muy alta.

  2. Iteración:

  • Se hace un cambio aleatorio a la solución actual \(S\) (buscar un “vecino”) \(S^*\).

  • Calcular el costo (o energía) de esa nueva solución \(S^*\).

  1. Regla de aceptación:

  • Si la nueva solución \(S^*\) es mejor, se acepta de inmediato.

  • Si la nueva solución \(S^*\) es peor, aun así podría aceptarse basándote en una probabilidad.

  1. Enfriamiento: conforme pasa el tiempo, la Temperatura (\(T\)) baja. Al bajar la temperatura, la probabilidad de aceptar soluciones peores disminuye.

  2. Estado Final: al final (cuando el metal está “frío”), el algoritmo deja de “saltar” y se asienta en la mejor solución que haya encontrado.

El buen funcionamiento reside en la probabilidad. Al principio, cuando la temperatura es alta, el algoritmo se mueve de forma casi aleatoria por todo el mapa de soluciones (exploración). Esto le permite saltar fuera de los “agujeros” o mínimos locales.

A medida que se enfría, se vuelve más selectivo y se dedica a perfeccionar la zona donde encontró los mejores resultados (explotación).

Resumen del funcionamiento

  1. Estado actual: solución \(S\) (rutas). Generación de solución inicial aleatoria o a partir de heurísticas.
  2. Vencindario \(S(N)\), soluciones cercanas, se pueden generar mediante búsqueda local.
  3. Temperatura \(T\): controla la probabilidad de aceptar soluciones peores.
  4. Regla de aceptación: si el cambio de costo es:

\[ \Delta = f(S^*) - f(S) \]

  • Si \(\Delta < 0\): se acepta (mejora)

  • Si \(\Delta >0\) se acepta con probabilidad:

\[ P = e^{\frac{-\Delta}{T}} \]

A medida que \(T\) baja existe menos probabilidad de aceptar empeoramientos.

Gráfico de Probabilidad de aceptación

Con Temperatura inicial grande:

  • Más exploración

  • Más probabilidad de escapar de óptimos locales

Con Temperatura inicial pequeña:

  • Comportamiento casi greedy desde el inicio

Ejemplo Recocido simulado para problema VRP en R

Durante las jornadas electorales en el Departamento de Córdoba, la Registraduría Nacional del Estado Civil debe realizar la recolección del material electoral sensible (urnas, formularios, actas y votos físicos) desde todos los municipios hacia el centro de consolidación ubicado en Montería.

Este proceso debe efectuarse inmediatamente después del cierre de las mesas de votación, debido a que:

  • El material es confidencial y de alto valor legal.

  • Debe llegar lo más rápido posible al centro de cómputo departamental.

  • Se busca minimizar riesgos de pérdida, manipulación o retrasos.

  • Los recursos de transporte y seguridad son limitados.

Para esta operación se dispone de una flota de \(12\) vehículos oficiales escoltados, cada uno con capacidad máxima de \(1.5\) toneladas, los cuales parten desde Montería, visitan un subconjunto de municipios para recolectar el material electoral y posteriormente regresan al centro de consolidación.

Cada municipio genera una cantidad estimada de material electoral (expresada en peso o número de cajas), por lo que la asignación de municipios a cada vehículo debe respetar la capacidad máxima de transporte. El peso del material electoral por municipio se encuentra disponible en el siguiente enlace: Peso de material electoral por municipio.

Adicionalmente, se dispone de una matriz de distancias en kilómetros entre municipios: Matriz de distancias en km.

Actualmente, las rutas se planifican de forma empírica, lo que genera:

  • mayores distancias recorridas,

  • incremento en los tiempos de consolidación de resultados,

  • mayores costos de combustible,

  • y mayores riesgos de seguridad por recorridos innecesarios.

En consecuencia, se requiere diseñar un modelo de ruteo de vehículos capacitado (Vehicle Routing Problem, VRP) que permita determinar:

  • qué municipios debe atender cada vehículo,

  • el orden óptimo de visita en cada ruta,

  • y minimizar la distancia total recorrida por la flota.

Datos

library(readxl)
library(dplyr)

#Se lee el archivo de Excel que contiene la matriz de distancias entre municipios
distancias_df<- read_excel("distancias_municipios_cordoba.xlsx")

#Matriz de distancias
municipios <- distancias_df[[1]]  #Extraer la primera columna que contiene los nombres de los municipios
D <- as.matrix(distancias_df[,-1]) #Construir la matriz de distancias eliminando la primera columna
rownames(D) <- municipios #Asignar nombres de filas y columnas a la matriz
colnames(D) <- municipios

#dataframe de demandas
demandas_df <- read_excel("demanda.xlsx")

#Matriz de demanda
demandas <- demandas_df$demanda #Crear vector de demandas
names(demandas) <- demandas_df$Municipio #Asignar nombres al vector de demandas para acceder por municipio

#Capacidad del vehículo (toneladas)
Q <- 1.5

#Nodo depósito (punto de origen y regreso de las rutas)
deposito <- "Monteria"

#Tamaño de flota (Número máximo de vehículos disponibles)
K <- 12

Solución inicial Aleatoria

#Define una función para generar rutas iniciales factibles para un VRP
generar_rutas_iniciales <- function(D, demandas, K, Q, deposito, max_intentos = 100) {
  
  #Obtiene los clientes excluyendo el depósito
  clientes_original <- setdiff(names(demandas), deposito)
  
  #Bucle que intenta generar una solución factible hasta max_intentos veces
  for (intento_global in 1:max_intentos) {
    
    #Genera una permutación aleatoria de los clientes
    clientes <- sample(clientes_original)
    
    #Crea una lista vacía de tamaño K (una ruta por vehículo)
    rutas <- vector("list", K)
    
    #Inicializa un vector con la carga actual de cada vehículo (empieza en 0)
    cargas <- rep(0, K)
    
    #Cada ruta empieza y termina en el depósito (sin clientes aún)
    for (k in 1:K) {
      rutas[[k]] <- c(deposito, deposito)
    }
    
    #Indicador de si la solución actual es factible
    factible <- TRUE
    
    #Itera sobre los clientes en orden aleatorio
    for (cliente in clientes) {
      
      #Genera un orden aleatorio de vehículos para intentar asignar el cliente
      vehiculos <- sample(1:K)
      
      #Indicador de si el cliente fue asignado a algún vehículo
      asignado <- FALSE
      
      #Intenta insertar el cliente en cada vehículo
      for (k in vehiculos) {
        
        #Verifica si al agregar el cliente no se excede la capacidad
        if (cargas[k] + demandas[cliente] <= Q) {
          
          #Inserta el cliente justo antes del último nodo (depósito final)
          rutas[[k]] <- append(rutas[[k]],
                               cliente,
                               after = length(rutas[[k]]) - 1)
          
          #Actualiza la carga del vehículo
          cargas[k] <- cargas[k] + demandas[cliente]
          
          #Marca que el cliente ya fue asignado
          asignado <- TRUE
          break
        }
      }
      
      #Si no se pudo asignar este cliente, romper intento
      if (!asignado) {
        
        #Marca la solución como no factible
        factible <- FALSE
        break
      }
    }
    
    #Si fue factible, devolver solución
    if (factible) {
      return(rutas)
    }
  }
  
  stop("No se pudo generar solución factible tras varios intentos")
}

Ejecutar función solución inicial

rutas <- generar_rutas_iniciales(D, demandas, K, Q, deposito)

Cantidad de rutas

length(rutas)

Cargas de cada ruta

sapply(rutas, function(r) sum(demandas[r]))

Costo de las rutas

distancia_ruta <- function(ruta, D){
  
  distancia <- 0
  
  for(i in 1:(length(ruta)-1)){
    
    distancia <- distancia + D[ruta[i], ruta[i+1]]
    
  }
  
  return(distancia)
}

distancias <- sapply(rutas, distancia_ruta, D)
distancia_total <- sum(distancias)
distancia_total

Resumen de resultados

data.frame(
  Ruta = 1:length(rutas),
  Distancia_km = sapply(rutas, distancia_ruta, D),
  Carga = sapply(rutas, function(r) sum(demandas[r]))
)

Función recocido simulado

#Función de Recocido Simulado para el VRP
simulated_annealing_vrp <- function(rutas, D, demandas, Q, T0 = 1000, alpha = 0.95, T_min = 1, iter_por_T = 100) {
  
  #Función costo
  costo_total <- function(rutas) {
    
    #Inicializa el costo total
    total <- 0
    
    #Recorre cada ruta
    for (r in rutas) {
      if (length(r) > 1) {
        for (i in 1:(length(r)-1)) {
          total <- total + D[r[i], r[i+1]]
        }
      }
    }
    return(total)
  }
  
  #Carga de ruta
  carga_ruta <- function(r) {
    sum(demandas[r])
  }
  
  #Ooperador: relocation
  vecino_relocation <- function(rutas) {
    
    #Copia de las rutas actuales
    rutas_new <- rutas
    
    #Selecciona dos rutas aleatorias
    r1 <- sample(1:length(rutas), 1)
    r2 <- sample(1:length(rutas), 1)
    
    #Si la ruta tiene solo depósito-cliente-depósito (o menos), no se puede extraer
    if (length(rutas[[r1]]) <= 3) return(rutas)
    
    #Selecciona un cliente (no depósito)
    i <- sample(2:(length(rutas[[r1]])-1), 1)
    
    #Cliente a mover
    cliente <- rutas[[r1]][i]
    
    #Elimina el cliente de la ruta origen
    rutas_new[[r1]] <- rutas[[r1]][-i]
    
    #Selecciona posición de inserción en la ruta destino
    pos <- sample(2:length(rutas[[r2]]), 1)
    
    #Inserta el cliente en la nueva ruta
    rutas_new[[r2]] <- append(rutas[[r2]], cliente, after = pos-1)
    
    #Si viola capacidad, se descarta el vecino
    if (carga_ruta(rutas_new[[r2]]) > Q) {
      return(rutas)
    }
    
    #Retorna la nueva solución vecina
    return(rutas_new)
  }
  
  # Operador swap
  vecino_swap <- function(rutas) {
    
    #Copia de rutas
    rutas_new <- rutas
    
    #Selecciona dos rutas aleatorias
    r1 <- sample(1:length(rutas), 1)
    r2 <- sample(1:length(rutas), 1)
    
    #Si alguna ruta no tiene clientes intercambiables, no se hace nada
    if (length(rutas[[r1]]) <= 3 || length(rutas[[r2]]) <= 3) {
      return(rutas)
    }
    
    #Selecciona posiciones de clientes en cada ruta
    i <- sample(2:(length(rutas[[r1]])-1), 1)
    j <- sample(2:(length(rutas[[r2]])-1), 1)
    
    #Clientes a intercambiar
    cliente1 <- rutas[[r1]][i]
    cliente2 <- rutas[[r2]][j]
    
    #Realiza el intercambio
    rutas_new[[r1]][i] <- cliente2
    rutas_new[[r2]][j] <- cliente1
    
    #Si alguna ruta viola capacidad, descarta el vecino
    if (carga_ruta(rutas_new[[r1]]) > Q ||
        carga_ruta(rutas_new[[r2]]) > Q) {
      return(rutas)
    }
    
    #Retorna solución vecina
    return(rutas_new)
  }
  
  #Inicialización simulated anealling
  
  #Temperatura inicial
  T <- T0
  
  #Solución actual
  actual <- rutas
  
  #Mejor solución encontrada
  mejor <- rutas
  
  #Costo de la solución actual
  costo_actual <- costo_total(actual)
  
  #Inicializa mejor costo
  costo_mejor <- costo_actual
  
  #Almacenamiento métrica
  historial_costo <- c() #Guarda costo de la solución actual en cada iteración
  historial_mejor <- c() #Guarda el mejor costo histórico
  historial_temp <- c() #Guarda la temperatura en cada iteración
  tasa_aceptacion <- c() #Guarda la tasa de aceptación por temperatura
  
  #Contador global de iteraciones
  iter_global <- 1
  
  #Ciclo principal
  
  #Ejecuta hasta que la temperatura sea suficientemente baja (enfriamiento)
  while (T > T_min) {
    
    #Contador de soluciones aceptadas en esta temperatura
    aceptados <- 0
    
    #Iteraciones internas a temperatura fija
    for (iter in 1:iter_por_T) {
      
      #Elegir operador
      if (runif(1) < 0.5) {
        vecino <- vecino_relocation(actual) #Con probabilidad 0.5 usa relocation
      } else {
        vecino <- vecino_swap(actual) #Con probabilidad 0.5 usa swap
      }
      
      costo_vecino <- costo_total(vecino) #Calcula costo del vecino
      delta <- costo_vecino - costo_actual #Diferencia de costos
      
      #Acepta si mejora o con probabilidad si empeora
      if (delta < 0 || runif(1) < exp(-delta / T)) {
        
        actual <- vecino #Actualiza solución actual
        costo_actual <- costo_vecino #Actualiza costo
        aceptados <- aceptados + 1 #Incrementa contador de aceptación
        
        #Actualiza "mejor" solución encontrada
        if (costo_actual < costo_mejor) {
          mejor <- actual
          costo_mejor <- costo_actual
        }
      }
      
      #Guardar historial global
      historial_costo[iter_global] <- costo_actual
      historial_mejor[iter_global] <- costo_mejor
      historial_temp[iter_global] <- T
      
      #Incrementa contador global
      iter_global <- iter_global + 1
    }
    
    #Tasa de aceptación por temperatura
    tasa_aceptacion <- c(tasa_aceptacion, aceptados / iter_por_T)
    
    #Enfriamiento, reduce la temperatura
    T <- alpha * T
  }
  
  #Retorna resultados del algoritmo
  return(list(
    mejor_rutas = mejor,
    mejor_costo = costo_mejor,
    historial_costo = historial_costo,
    historial_mejor = historial_mejor,
    historial_temp = historial_temp,
    tasa_aceptacion = tasa_aceptacion
  ))
}

Ejecución del algoritmo

resultados <- simulated_annealing_vrp(rutas, D, demandas, Q, T0=1000, alpha=0.95, T_min=1)

Gráficos

plot(resultados$historial_costo, type = "l",
     main = "Evolución del costo",
     xlab = "Iteración", ylab = "Costo")

lines(resultados$historial_mejor, col = "red", lwd = 2)
legend("topright", legend = c("Actual", "Mejor"),
       col = c("black", "red"), lty = 1)

Cantidad de rutas

length(resultados$mejor_rutas)

Cargas de cada ruta

sapply(resultados$mejor_rutas, function(r) sum(demandas[r]))

Costo de las rutas

costo_sa <- resultados$mejor_costo

Resumen de resultados

data.frame(
  Ruta = 1:length(resultados$mejor_rutas),
  Distancia_km = sapply(resultados$mejor_rutas, distancia_ruta, D),
  Carga = sapply(resultados$mejor_rutas, function(r) sum(demandas[r]))
)