Problemas combinatorios, heurísticas y metaheurísticas

Problemas de optimización.

En el contexto científico la optimización es el proceso de tratar de encontrar la mejor solución posible para un determinado problema. En un problema de optimización existen diferentes soluciones, un criterio para discriminar entre ellas y el objetivo es encontrar la mejor. De forma más precisa, estos problemas se pueden expresar como encontrar el valor de unas variables de decisión para los que una determinada función objetivo alcanza su valor máximo o mínimo. El valor de las variables en ocasiones está sujeto a unas restricciones.

Se pueden encontrar una gran cantidad de problemas de optimización estudiados clásicamente en ingeniería industrial y aplicados a problemas o decisiones logísticas (servicio al cliente o procesamiento de pedidos, inventarios, transporte y localización).

1. Problemas de asignación de recursos (AP = Assignment Problem).

Problema clásico de optimización cuyo propósito es hallar la forma de asignar un conjunto \(I\) de recursos disponibles (personas, máquinas, robots, computadores, instalaciones…) para la realización de un conjunto \(J\) de actividades o de lugares (puestos de trabajo, operaciones productivas, tareas de ensamble, tareas computacionales… ) al menor costo posible y suponiendo que cada recurso se destina a una única actividad o lugar y que cada actividad o lugar solo puede ser ejecutada por uno de los recursos. Se supone, por tanto que el número de recursos disponibles es igual al número de actividades a realizar o de lugares a ocupar, siendo tal número igual a \(n= |I| = |J|\).

Por ejemplo, suponga un AP con \(n=4\), como se observa en la siguiente figura:

El número de asignaciones posibles es \(4!=24\).

El AP en cuestión, a pesar de la simpleza de su enunciado, tiene numerosas aplicaciones prácticas dentro de la Ingeniería Industrial, y específicamente dentro de las decisiones logísticas, ya que lo de tener que decidir cómo repartir unidades de recursos productivos entre unidades de actividades productivas, cuando los costos de asignación son lineales o pueden suponerse como tales, forma parte de muchas decisiones habituales en la prácticas. Entre las aplicaciones del problema de asignación se encuentran las siguientes:

• Asignar instalaciones a localizaciones.

• Asignar cuotas de producción de fábricas a centros de distribución con demanda.

• Asignar stocks de productos a almacenes.

• Emparejar vehículos con rutas de reparto o de recogida.

• Definir dónde instalar bodegas para e-commerce

• Asignar centros de distribución a regiones.

El AP se puede formular matemáticamente como sigue:

• Conjuntos

  • \(I\): conjunto de recursos, \(I=\{1,2,...,n \}\)
  • \(J\): conjunto de actividades o lugares, \(J=\{1,2,...,n \}\)

Se supone que: \(|I| = |J| = n\)

• Parámetros

\(c_{ij}\): costo de asignar el recurso \(i \in I\) a la actividad o lugar \(j \in J\).

• Variables de decisión

\[\begin{align} x_{ij} = \begin{cases} 1:~ \text{si el recurso }i \in I \text{ se asigna a actividad o lugar } j \in J\\ 0:~ \text{dlc } \end{cases} \end{align}\]

• Función objetivo

\[ min Z = \sum_{i \in I} \sum_{j \in J} c_{ij}x_{ij} \]

• Restricciones

Cada recurso se asigna a una única actividad o lugar:

\[ \sum_{j \in J} x_{ij} = 1~~~ \forall i \in I \]

Cada actividad o lugar recibe un único recurso:

\[ \sum_{i \in I} x_{ij} = 1~~~ \forall j \in J \]

Naturaleza de la variable

\[ x_{ij} \in \{0,1 \} ~~~ \forall i \in I, j \in J \]

Ejemplo problema AP en R

Una empresa planea ubicar 4 fábricas en 4 parques industriales disponibles. Cada combinación fábrica–parque tiene un costo logístico anual estimado, que incluye:

• transporte a mercados,

• acceso a proveedores,

• costos operativos del parque industrial.

Los costos asociados (en miles de millones de pesos/año) son los siguientes:

Fábrica/Parque	P1	P2	P3	P4
F1	14	5	8	7
F2	2	12	6	5
F3	7	8	3	9
F4	2	4	6	10

El modelo matemático sería

Conjuntos

• \(I= \{1,2,3,4 \}\): fábricas
• \(J= \{1,2,3,4 \}\): parques industriales

Parámetro

\(c_{ij}\): costo de ubicar fábrica \(i \in I\) en el parque industrial \(j \in J\) (Tabla de costos).

Variable de decisión

\[\begin{align} x_{ij} = \begin{cases} 1:~ \text{si la fábrica } i ~\text{se asigna al parque industrial } j \\ 0:~ \text{dlc } \end{cases} \end{align}\]

• Función objetivo

\[ min~Z = \sum_{1}^4 \sum_{1}^4 c_{ij}x_{ij} \]

• Restricciones

Cada fábrica se asigna a un solo parque industrial:

\[ \sum_{1}^4 x_{ij} = 1~~~ \forall i \in I \]

Cada parque industrial recibe una única fábrica:

\[ \sum_{1}^4 x_{ij} = 1~~~ \forall j \in J \]

Naturaleza de la variable

\[ x_{ij} \in \{0,1 \} ~~~ \forall i \in I, j \in J \]

Librerías necesarias

# Cargar librerías
library(ompr) #Formular modelos LP, MIP
library(ompr.roi) #Conecta modelo con solver
library(ROI) #Infraestructura para solver
library(ROI.plugin.symphony) #Solver
library(magrittr) #Escritura del modelo

Definición de datos

# Parámetro costos
costos <- matrix(c(14, 5, 8, 7, 2, 12, 6, 5, 7, 8, 3, 9, 2, 4, 6, 10),
  nrow = 4,
  byrow = TRUE
)

n <- nrow(costos)

Modelo

# Modelo
modelo <- MIPModel() %>%
  
  # Variable binaria x[i,j]
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
  
  # Función objetivo
  set_objective(
    sum_expr(costos[i, j] * x[i, j], i = 1:n, j = 1:n),
    "min"
  ) %>%
  
  # Cada fábrica a un solo parque industrial
  add_constraint(
    sum_expr(x[i, j], j = 1:n) == 1,
    i = 1:n
  ) %>%
  
  # Cada parque industrial recibe una sola fábrica
  add_constraint(
    sum_expr(x[i, j], i = 1:n) == 1,
    j = 1:n
  )

Resolver con solver

# resultado 
resultado <- solve_model(
  modelo,
  with_ROI(
    solver = "symphony",
    verbosity = 1
  )
)

Extracción de la solución

# extracción de la solución
solucion <- get_solution(resultado, x[i, j]) %>%
  dplyr::filter(value > 0.5)
solucion

2. Problema de transporte (TP: transportation problem).

Una evidente generalización del AP consiste en suponer que tanto los recursos como las actividades o lugares están replicados (Un recurso puede tener capacidad mayor que 1, una actividad puede requerir varias unidades de recurso).

Suponga que el tipo de recurso \(i \in I\) corresponde a una fábrica de productos que tiene asociada una oferta de unidades igual a \(o_i~\forall i\) y que todo tipo de actividad o lugar \(j \in J\) corresponde a un centro de consumo de tales productos, presentando una demanda \(d_j~\forall j\).

Suponga, además, que el costo de llevar una unidad de producto desde la fábrica \(i \in I\) hasta el centro de consumo \(j \in J\) es \(c_{ij}~ \forall i, \forall j\). Ante este escenario el número de tipos de recursos o fábricas (\(n= |i|\)) no siempre coincide con el número de tipos de actividad o centros de consumo (\(m= |j|\)).

Por su parte, ahora la relación entre la oferta de las fábricas y las demandas de los centros de consumo puede responder a las tres situaciones siguientes:

1. La suma de ofertas coincide con la suma de demandas, esto es, \(\sum_{i \in I} o_i = \sum_{j \in J} d_j\), en este caso el problema se considera equilibrado.

2. La suma de ofertas es mayor a la suma de demandas, esto es, \(\sum_{i \in I} o_i > \sum_{j \in J} d_j\). En este caso se puede definir y añadir al conjunto \(J\) un centro ficticio de consumo \(\omega\), al que se le asigna una demanda ficticia igual al exceso de oferta, por lo tanto:

\[d_{\omega} = \sum_{i \in I} o_i - \sum_{j \in J} d_j\]

Además, se le pueden agregar unos costos ficticios hasta \(\omega\) iguales a cero, por lo tanto:

\[c_{i \omega} = 0,~ \forall i\]

Con lo anterior, el problema queda equilibrado.

3. La suma de ofertas es menor a la suma de demandas, esto es, \(\sum_{i \in I} o_i < \sum_{j \in J} d_j\). Para lo anterior se definirá una fábrica ficticia \(\alpha\) igual a la brecha de demanda:

\[ o_{\alpha} = \sum_{j \in J} d_j - \sum_{i \in I} o_i \]

Además, se le pueden agregar unos costos de transporte ficticios desde \(\alpha\) iguales a cero, por lo tanto:

\[c_{\alpha j} = 0,~ \forall j\]

Con lo anterior, el problema queda equilibrado.

La descripción anterior, corresponde a un problema de asignación de recursos conocido como Problema de Transporte (TP: transportation problem). Una formulación matemática posible es la siguiente:

• Conjuntos

  • \(I\): conjunto de fábricas, \(I= \{ 1,2,...,n \}\)
  • \(J\): conjunto de centros de consumo, \(J= \{1,2,...,m \}\)

• Parámetros

  • \(o_i\): oferta disponible en la fábrica \(i \in I\) (Enteras).
  • \(d_j\): demanda del centro de consumo \(j \in J\) (Enteras).
  • \(c_{ij}\): costo de transportar una unidad desde \(i \in I\) hasta \(j \in J\).

• Variable de decisión

\[ x_{ij} \ge 0~\forall i \in I, \forall j \in J \\ x_{ij} \in \mathbb{Z} \]

Cantidad de producto transportado desde la fábrica \(i \in I\) hasta el centro de consumo \(j \in J\).

• Función objetivo

\[ min~Z = \sum_{i \in I} \sum_{j \in J} c_{ij}x_{ij} \]

• Restricciones

De oferta: la oferta de cada fábrica debe ser distribuida completamente:

\[ \sum_{j \in J} x_{ij} = o_{i}, ~~~\forall~i \in I \]

De demanda: cada centro de consumo debe recibir su demanda.

\[ \sum_{i \in I} x_{ij} = d_{j}, ~~~\forall~j \in J \]

Naturaleza de la variable

\[x_{ij} \ge 0~~~ \forall i \in I, \forall j \in J \\ x_{ij} \in \mathbb{Z}\]

La restricción \(x_{ij} \in \mathbb{Z}\) no es estrictamente necesaria.

El AP es un caso particular del TP, basta con imponen \(o_i = d_j = 1~\forall i \in I, ~\forall j \in J\)

Ejemplo problema TP en R

Una empresa productora de cemento tiene 3 plantas de producción ubicadas en distintas regiones del país y debe abastecer 4 centros de distribución regionales.

Cada planta tiene una capacidad limitada de producción y cada centro de distribución tiene una demanda mínima que debe ser satisfecha completamente.

El costo de transporte depende de la distancia, peajes y tipo de vía entre cada planta y cada centro. El objetivo de la empresa es minimizar el costo total de transporte, cumpliendo todas las demandas sin exceder la oferta disponible

La oferta por planta (en miles de toneladas) se muestra en la siguiente tabla:

Planta	Oferta
F1	80
F2	120
F3	100

La demanda de los centros de distribución (en miles de toneladas) son las siguientes:

Centro Distribución	Demanda
C1	70
C2	90
C3	80
C4	60

Los costos del transporte está dado en la siguiente tabla:

Planta	Centro 1	Centro 2	Centro 3	Centro 4
F1	4	6	9	5
F2	5	4	7	6
F3	6	7	4	3

Conjuntos

• \(I= \{1,2,3 \}\): plantas de produccción de cemento.
• \(J= \{1,2,3,4 \}\): centros de distribución.

Parámetros (Tablas)

• \(o_i\): oferta disponible de la planta productora \(i \in I\).
• \(d_j\): demanda del centro de distribución \(j \in J\).
• \(c_{ij}\): costo de transportar desde la planta productora \(i \in I\) hasta el centro de distribución \(j \in J\).

Variable de decisión

\[ x_{ij} \ge 0 \]

Cantidad de cemento transportado desde la planta \(i \in I\) al centro de distribución \(j \in J\).

Función objetivo

\[ min~Z = \sum_{i=1}^3 \sum_{j=1}^4 c_{ij}x_{ij} \]

Restricciones

Oferta

\[ \sum_{j=1}^4 x_{ij} = o_i,~\forall i \]

Demanda

\[ \sum_{i=1}^3 x_{ij} = d_j,~\forall j \]

Naturaleza de la variable

\[x_{ij} \ge 0~~~ \forall i \in I, \forall j \in J \]

Librerías necesarias

# Cargar librerías
library(ompr)
library(ompr.roi)
library(ROI)
library(ROI.plugin.symphony)
library(magrittr)

Definición de datos

costos <- matrix(c(4, 6, 9, 5, 5, 4, 7, 6, 6, 7, 4, 3), nrow = 3, byrow = TRUE )

# Oferta y demanda
oferta  <- c(80, 120, 100)
demanda <- c(70, 90, 80, 60)

n <- length(oferta)   # plantas
m <- length(demanda)  # centros

Modelo

modelo_tp <- MIPModel() %>%
  
  # Variable continua x[i,j] ≥ 0
  add_variable(x[i, j], i = 1:n, j = 1:m, lb = 0) %>%
  
  # Función objetivo
  set_objective(
    sum_expr(costos[i, j] * x[i, j], i = 1:n, j = 1:m),
    "min"
  ) %>%
  
  # Restricciones de oferta
  add_constraint(
    sum_expr(x[i, j], j = 1:m) == oferta[i],
    i = 1:n
  ) %>%
  
  # Restricciones de demanda
  add_constraint(
    sum_expr(x[i, j], i = 1:n) == demanda[j],
    j = 1:m
  )

Resolver con solver

resultado_tp <- solve_model(
  modelo_tp,
  with_ROI(
    solver = "symphony",
    verbosity = 1
  )
)

Extracción de la solución

solucion_tp <- get_solution(resultado_tp, x[i, j]) %>%
  dplyr::filter(value > 0)

solucion_tp

3. Problema de asignación cuadrática (QAP: Quadratic Assignment Problem).

Suponga que se dispone de \(n\) instalaciones que realizan la misma actividad industrial, como por ejemplo fábricas, almacenes, estaciones de generación de energía, centros de distribución o torres de comunicación. Estas instalaciones conforman el conjunto \(I\).

Entre cada par de instalaciones existe algún tipo de interacción, la cual puede representarse mediante un flujo, como el flujo de productos, de energía o de información.

Además, se dispone de \(n\) posibles emplazamientos o ubicaciones, que conforman el conjunto \(J\). Entre estos emplazamientos existen distancias que no son despreciables y que influyen directamente en los costos asociados a la ubicación de las instalaciones.

En tal caso, los costos derivados de las comunicaciones entre instalaciones dependen, al menos, de los flujos \(\phi_{ik}\) entre parejas de instalaciones \(\{i,k\} \subseteq I\) y las distancias \(\delta_{jl}\) entre pares de emplazamientos \(\{ j,l \} \subseteq J\).

Los costos de comunicación, que se denotan \(c_{ijkl}\), se puede representar mediante una función \(f( \cdot )\) que puede ser no lineal: \(c_{ijkl} = f(\phi_{ik,~\delta_{jl}});~\forall \{i,k \} \subseteq I, ~\forall \{j,l \} \subseteq J\).

En tales condiciones, se está ante un problema clásico de localización denominado Problema de Asignación Cuadrática (QAP: Quadratic Assignment Problem). Consiste en hallar una asignación de las instalaciones del conjunto \(I\) a los emplazamientos del conjunto \(J\) al menor costo posible, suponiendo que cada instalación se localiza en un único emplazamiento y que cada emplazamiento solo puede ser ocupado por una instalación.

Una formulación matemática posible es la siguiente:

Conjuntos

• \(I\): conjunto de instalaciones, \(I=\{1,2,...,n\}\).
• \(J\): conjunto de emplazamientos o lugares, \(J= \{1,2,...,n \}\)

\(n\) representa el número total de instalaciones y emplazamientos.

Parámetro

\(c_{ijkl}\): costo asociado a la asignación de instalaciones \(i \in I\) y \(k \in I\) a los emplazamientos \(j \in J\) y \(l \in J\) respectivamente.

Variables de decisión

\[\begin{align} x_{ij} = \begin{cases} 1:~ \text{si la instalación }i \in I \text{ se asigna al emplazamiento o lugar } j \in J\\ 0:~ \text{dlc } \end{cases} \end{align}\]

Función objetivo

\[ min~Z = \sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n \sum_{l=1}^n c_{ijkl} \times x_{ij} \times x_{kl} \]

La función objetivo representa la minimización de los costos de comunicación entre instalaciones una vez han sido localizadas en el conjunto de emplazamientos.

Restricciones

Toda instalación \(i \in I\) debe ser asignada a un emplazamiento \(j \in J\)

\[ \sum_{j=1}^n x_{ij} = 1~~~\forall i=1,2,...,n \]

Todo emplazamiento \(j \in J\) debe ser ocupado por una instalación \(I \in I\)

\[ \sum_{i=1}^n x_{ij} = 1~~~\forall j=1,2,...,n \]

Naturaleza de la variable

\[\begin{align} x_{ij} \in \{ 0, 1\}~~~\forall i=1,2,...,n;~\forall j=1,2,...,n \end{align}\]

Ejemplo problema QAP en R

Por adquisión de predios una empresa de distribución regional que opera cuatro centros operativos que clasifican y redistribuyen mercancía debe decidir cómo ubicar estos cuatro centros en cuatro predios disponibles dentro de una zona geográfica determinada.

Entre los centros operativos se intercambia inventario para cumplir con los requerimientos de los clientes cuando sea necesario. El flujo \(\phi_{ik}\) que se ha dado entre ellos se muestra en la siguiente tabla (cifras en toneladas):

Centro \(i\) - Centro \(k\)	1	2	3	4
1	0	2	15	1
2	2	0	3	10
3	15	3	0	4
4	1	10	4	0

Por otra parte, los cuatro predios disponibles se encuentran separados por distancias geográficas que influyen directamente en los costos de transporte y comunicación. Las distancias (en cientos de kilómetros) se muestran en la siguiente tabla:

Localización \(j\) - Localización \(l\)	1	2	3	4
1	0	8	6	1
2	8	0	4	7
3	6	4	0	5
4	1	7	5	0

El costo del intercambio entre dos centros operativos depende simultáneamente:

• del flujo existente entre dichos centros, y

• de la distancia entre los emplazamientos en los que son ubicados.

En este caso, se asume que el costo asociado a asignar los centros \(i\) y \(k\) a los emplazamientos \(j\) y \(l\) está dado por:

\[ c_{ijkl} = \phi_{ik} \times \delta_{jl} \]

El objetivo de la empresa es determinar la mejor asignación de centros a localizaciones que minimice el costo total de interacción, bajo la condición de que:

• cada centro operativo se ubique en un único emplazamiento, y

• cada emplazamiento sea ocupado por un solo centro operativo.

Conjuntos

• \(I\): conjunto de centros operativos, \(I=\{1,2,3,4\}\).
• \(J\): conjunto de elocalizaciones, \(J= \{1,2,3,4 \}\)

Parámetro

• \(\phi_{ik}\): flujo entre los centros operativos \(i \in I\) y \(k \in I\)
• \(\delta_{ik}\):distancia entre las localizaciones \(j \in J\) y \(l \in J\)
• \(c_{ijkl}\): costo de ubicar centros operativos \(i \in I\) y \(k \in I\) a las localizaciones \(j \in J\) y \(l \in J\) respectivamente. \(c_{ijkl} = \phi_{ik} \times \delta_{jl}\).

Variables de decisión

\[\begin{align} x_{ij} = \begin{cases} 1:~ \text{si el centro operativo }i \in I \text{ se ubica en la localización } j \in J\\ 0:~ \text{dlc } \end{cases} \end{align}\]

Función objetivo

\[ min~Z = \sum_{i=1}^4 \sum_{j=1}^4 \sum_{k=1}^4 \sum_{l=1}^4 c_{ijkl} \times x_{ij} \times x_{kl} \]

Restricciones

Asignación única de centros

\[\sum_{j=1}^n x_{ij}=1;~~~\forall i \in I\]

Ocupación única de localizaciones

\[\sum_{i=1}^n x_{ij}=1;~~~\forall j \in J\]

Naturaleza de las variables

\[x_{ij} \in \{0,1\}\]

Librerías

# Librerías
library(ompr)
library(ompr.roi)
library(ROI)
library(ROI.plugin.symphony)
library(magrittr)

Datos del problema

# Matriz de flujos
phi <- matrix(c(0,  2, 15,  1, 2,  0,  3, 10, 15,  3,  0,  4, 1, 10,  4,  0), nrow = 4, byrow = TRUE)

# Matriz de distancias
delta <- matrix(c(0, 8, 6, 1, 8, 0, 4, 7, 6, 4, 0, 5, 1, 7, 5, 0), nrow = 4, byrow = TRUE)

n <- nrow(phi)

Construcción parámetro \(c_{ijkl}\)

cijkl <- array(0, dim = c(n, n, n, n))

for (i in 1:n)
  for (k in 1:n)
    for (j in 1:n)
      for (l in 1:n)
        cijkl[i, j, k, l] <- phi[i, k] * delta[j, l]

Modelo

modelo_qap <- MIPModel() %>%
    add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
    set_objective(sum_expr(cijkl[i, j, k, l] * x[i, j] * x[k, l],i = 1:n, j = 1:n, k = 1:n, l = 1:n),sense = "min") %>%
    add_constraint(sum_expr(x[i, j], j = 1:n) == 1,i = 1:n) %>%
    add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 1:n)

Resultado

resultado <- solve_model(modelo_qap, with_ROI(solver = "symphony"))

Extracción

solucion <- get_solution(resultado, x[i, j]) %>%
  dplyr::filter(value > 0.5)

solucion

Con la herramienta actual no se puede resolver, pues el problema no es lineal (\(c_{ijkl} \times x_{i, j} \times x_{k,l}\)) por lo que es necesario Linealizar.

Linealización

Se define nueva variable binaria:

\[ y_{ijkl} = x_{ij} \times x_{kl} \]

Se reemplaza el término cuadrático por \(y_{ijkl}\) y se agregan restricciones que hagan a la linealización equivalente al problema original.

Variables

• \(x_{ij} \in \{0,1\}\)
• \(y_{ijkl} \in \{0,1\}\)

Función objetivo lineal

\[ min~Z = \sum_{i,j,k,l} c_{ijkl} \times y_{ijkl} \]

Restricciones de linealización

\[\begin{align} y_{ijkl} &\le x_{ij} \\ y_{ijkl} &\le x_{kl} \\ y_{ijkl} &\ge x_{ij} + x_{kl} - 1\\ y_{ijkl} &\ge 0 \end{align}\]

Implementación con linealización

#Librerías
library(ompr)
library(ompr.roi)
library(ROI)
library(ROI.plugin.symphony)
library(magrittr)

#Datos
# Matriz de flujos
phi <- matrix(c(0,  2, 15,  1, 2,  0,  3, 10, 15,  3,  0,  4, 1, 10,  4,  0), nrow = 4, byrow = TRUE)

# Matriz de distancias
delta <- matrix(c(0, 8, 6, 1, 8, 0, 4, 7, 6, 4, 0, 5, 1, 7, 5, 0), nrow = 4, byrow = TRUE)

n <- nrow(phi)

# Construcción costo
cijkl <- array(0, dim = c(n, n, n, n))
for (i in 1:n)
  for (k in 1:n)
    for (j in 1:n)
      for (l in 1:n)
        cijkl[i, j, k, l] <- phi[i, k] * delta[j, l]

#Modelo linealizado
modelo_qap_lin <- MIPModel() %>%
  
  # Variables x
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
  
  # Variables y
  add_variable(y[i, j, k, l], 
               i = 1:n, j = 1:n, k = 1:n, l = 1:n,
               type = "binary") %>%
  
  # Función objetivo lineal
  set_objective(
    sum_expr(cijkl[i, j, k, l] * y[i, j, k, l],
             i = 1:n, j = 1:n, k = 1:n, l = 1:n),
    sense = "min"
  ) %>%
  
  # Asignación uno a uno
  add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n) %>%
  add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 1:n) %>%
  
  # Restricciones de linealización
  add_constraint(y[i, j, k, l] <= x[i, j],
                 i = 1:n, j = 1:n, k = 1:n, l = 1:n) %>%
  
  add_constraint(y[i, j, k, l] <= x[k, l],
                 i = 1:n, j = 1:n, k = 1:n, l = 1:n) %>%
  
  add_constraint(y[i, j, k, l] >= x[i, j] + x[k, l] - 1,
                 i = 1:n, j = 1:n, k = 1:n, l = 1:n)

#Resolver
resultado <- solve_model(
  modelo_qap_lin,
  with_ROI(solver = "symphony")
)

#Solución
solucion <-get_solution(resultado, x[i, j]) %>%
  dplyr::filter(value > 0.5)
solucion

• El modelo crece como \(O(n^4)\).
• Para \(n>8\) el enfoque resuelto no es práctico.
• Se resuelve con heurística o metaheurísticas.

4. Problema del agente viajero (TSP: traveling salesman problem).

Es uno de los problemas mayormente estudiados en investigación de operaciones. Cuando la Teoría de la Complejidad Algorítmica se desarrolló, el TSP fue uno de los primeros problemas en estudiarse probándose que pertenece a la clase NP-Hard. Desde los métodos de Ramificación y Acotamiento hasta los basados en la Combinatoria Poliédrica, pasando por los procedimientos Metaheurísticos, todos han sido inicialmente probados en el TSP, convirtiéndose ´este en una prueba obligada para validar cualquier técnica de resolución de problemas enteros o combinatorio.

El problema puede enunciarse como: un agente viajero de comercio debe visitar \(n\) ciudades, comenzando y finalizando en la misma ciudad. Conociendo el costo de de ir de una ciudad a otra, determinar el recorrido de costo mínimo.

Para enunciar el problema de manera formal, se define un grafo \(G(V,A,C)\), donde:

• \(V\): es el conjunto de vértices.
• \(A\): es el conjunto de aristas.
• \(C\): \(c_{ij}\) es el costo o distancia de la arista \((i,j)\).

Además, se define:

• Un camino es una sucesión de aristas, en donde el vértice final de cada arista coincide con el final de la siguiente. Se puede representar por la sucesión de vértices utilizados.

• Un camino es simple si no utiliza el mismo vértice más de una vez.

• Un ciclo es una camino en donde el vértice final coincide con el inicial.

• Un ciclo es simple si lo es el camino que lo define.

• Un subtour es un ciclo simple que no pasa por todo los vértices del grafo.

• Un tour o ciclo hamiltoniano es un ciclo simple que pasa por todos los vértices.

El TSP consiste en determinar un tour o ciclo hamiltoniano a mínimo costo. Se suele utilizar para:

• Ruteo de vehículos.
• Recogida de material en almacenes (Robotizado).
• Recolección de paquetería.
• Picking dentro de almacenes.
• Programación de visitas técnicas.

Un planteamiento matemático es el siguiente:

Conjuntos

• \(V\): conjunto de vértices (ciudades, nodos, lugares, localizaciones). \(V=\{1,2,...,n \}\)

Parámetros

• \(c_{ij}\): costo de ir desde el vértice \(i \in V\) al vértice \(j \in V\) (puede ser heterogéneo).

Variables de decisión

\[ x_{ij} \in \{0,1 \} \\ x_{ij} = \begin{cases} 1, & \text{si se viaja de } i \in V \text{ a } j \in V \\ 0, & \text{dlc } \end{cases} \]

Restricciones

• De cada vértice \(i \in V\) sale exactamente un arco.

\[ \sum_{j \in V, ~j \ne i} x_{ij} =1; \forall i \in V \]

• A cada vértice \(i \in V\) entra exactamente un arco,

\[ \sum_{i \in V, ~i \ne j} x_{ij} =1; \forall j \in V \]

• Eliminación de subtours (MTZ: Miller-Tucker-Zemlin)

Variable extra: \(u_i\): orden de visita.

\(u_i\): posición del vértice \(i \in V\) en la ruta. \(1\leq u_i \leq n\) (aplica para \(i \ge 2\)).

Entonces, para \(i \ne j\) y \(i,j \in \{2,...,n \}\)

\[ u_i - u_j + nx_{ij} \leq n-1 \]

En conclusión:

\[\begin{align} &u_i - u_j + nx_{ij} \leq n-1;~ \forall i \ne j, i,j \ge 2\\ &1\leq u_i \leq n \\ \end{align}\]

Naturaleza de las variables

\[\begin{align} x_{ij} &\in \{0,1 \} \\ &u_i \ge 0 \end{align}\]

Ejemplo problema TSP 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 después del cierre de las mesas de votación, ya 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 son limitados.

Para este proceso se dispone de un vehículo oficial escoltado, el cual parte desde Montería, debe visitar cada municipio para recoger el material y regresar a Montería.

Actualmente, las rutas se planifican de forma empírica, lo que genera mayores distancias recorridas,aumento en tiempos de consolidación de resultados, mayores costos de combustible y riesgos de seguridad por recorridos innecesarios.

Para lo anterior se dispone de una matriz de distancias entre municipios que se peude ver en el siguiente enlace: Matriz de distancias en km.

Conjuntos

• $V: {1,2,3,…,30 } $: conjunto de municipios.
  • Nodo 1: Montería (origen y destino).

Parámetros

• \(c_{ij}\): distancia por carretera desde el municipio \(i \in V\) hasta el municipio \(j \in V\) (\(c_{ij} \ge 0,~c_{ii} = 0\))

Variables de decisión

• Variable binaria de la ruta

\[ x_{ij} = \begin{cases} 1, & \text{si el vehículo viaja de } i \in V \text{ a } j \in V \\ 0, & \text{dlc } \end{cases} \]

• Variable auxiliar MTZ, posición del municipio \(i \in V\) en la secuencia de visita.

\[\begin{align} u_i \in \mathbb{R} \end{align}\]

Función objetivo

\[\begin{align} min~Z: \sum_{i \in V} \sum_{j \in V} c_{ij}x_{ij} \end{align}\]

Restricciones

• Salir exáctamente una vez de cada municipio

\[\begin{align} \sum_{j \in V,~i \ne j} x_{ij} = 1,~ \forall i \in V \end{align}\]

• Llegar exactamente una vez a cada municipio

\[\begin{align} \sum_{i \in V,~i \ne j} x_{ij} = 1,~ \forall j \in V \end{align}\]

• Eliminación de subtours (MTZ)

Para todos \(i \ne j;~~i,j \in \{2,3,4,...,30 \}\)

\[\begin{align} u_i - u_j +30x_{ij} \leq 29 \end{align}\]

• Fijación de origen en Montería

\[\begin{align} u_i = 1 \end{align}\]

• Cotas de orden

\[\begin{align} 2 \leq u_i \leq 30;~ \forall i=2,3,...,30 \end{align}\]

• Naturaleza de las variables

\[\begin{align} x_{ij} \in \{0,1 \} \\ u_i \ge 0 \end{align}\]

Librerías necesarias

# Cargar librerías
library(ompr)
library(ompr.roi)
library(ROI)
library(ROI.plugin.symphony)
library(readxl)
library(dplyr)
library(leaflet)

Definición de datos

#Dataframe de distancias
distancias_df <- read_excel("distancias_municipios_cordoba.xlsx")

  #Nombres de municipios
municipios_vec <- distancias_df[[1]]

  #Matriz de distancias
dist_matrix <- as.matrix(distancias_df[, -1])

  #Asignación de nombres a matriz de distancias
rownames(dist_matrix) <- municipios_vec
colnames(dist_matrix) <- municipios_vec

#Tamaño del problema
n <- nrow(dist_matrix)

#Data frame con coordenadas
municipios_df <- read_excel("C:/Users/000322041/OneDrive - UPB/UPB/Asignaturas_2026_1/Metodos_cuantitativos_logistica_27043/Clases/0_Combinatorios_heuristicas_metaheuristicas/coordenadas.xlsx", 
                            sheet = "Hoja2")
  #orden debe municipios df debe coincidir con municipios_vec

#Nodo inicial
nodo_inicial <- which(municipios_vec == "Monteria")

Modelo

modelo <- MIPModel() %>%
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
  add_variable(u[i], i = 1:n, lb = 1, ub = n) %>%
  set_objective(sum_expr(dist_matrix[i, j] * x[i, j], i = 1:n, j = 1:n), "min") %>%
  add_constraint(x[i, i] == 0, i = 1:n) %>% #No viajar de un nodo a sí mismo
  add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = 1:n) %>%   #Una salida por nodo
  add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = 1:n) %>%   # Una entrada por nodo
  add_constraint(u[nodo_inicial] == 1) %>% # Fijar Montería como nodo inicial
  add_constraint(u[i] - u[j] + n * x[i, j] <= n - 1, i = 2:n, j = 2:n, i != j)

Resolver el modelo

resultado <- solve_model(
  modelo,
  with_ROI(
    solver = "symphony",
    control = list(time_limit = 3600),
    verbosity = -1
  )
)

Extracción de la solución

#Solución
solucion <- get_solution(resultado, x[i, j]) %>%
  filter(value > 0.5)

#Ruta
ruta_idx <- integer(n)
ruta_idx[1] <- nodo_inicial

for (k in 2:n) {
  ruta_idx[k] <- solution$j[solution$i == ruta_idx[k - 1]]
}

ruta_idx <- c(ruta_idx, ruta_idx[1])

#Ruta con nombres
ruta_municipios <- municipios_vec[ruta_idx]
ruta_municipios

#Costo total
dist_total <- sum(
  dist_matrix[
    cbind(ruta_idx[-length(ruta_idx)], ruta_idx[-1])
  ]
)

Mapa de ruta

#Dataframe ordenado para el mapa
ruta_df <- municipios_df[ruta_idx, ]
ruta_df$orden <- 1:nrow(ruta_df)

#Mapa básico
leaflet(ruta_df) %>%
  addTiles() %>%
  addMarkers(
    lng = ~longitud,
    lat = ~latitud,
    popup = ~Municipio
  ) %>%
  addPolylines(
    lng = ~longitud,
    lat = ~latitud,
    weight = 4
  )

#Mapa con orden

leaflet(ruta_df) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~longitud,
    lat = ~latitud,
    label = ~paste0(orden, ". ", Municipio),
    radius = 6
  ) %>%
  addPolylines(
    lng = ~longitud,
    lat = ~latitud
  )

Código creación matriz distancias

#Librerías
library(httr)
library(jsonlite)
library(dplyr)
library(openxlsx)


#Datos de coordenadas geográficas
library(readxl)
municipios <- read_excel("C:/Users/000322041/OneDrive - UPB/UPB/Asignaturas_2026_1/Metodos_cuantitativos_logistica_27043/Clases/0_Combinatorios_heuristicas_metaheuristicas/coordenadas.xlsx", 
                          sheet = "Hoja2")

#API Key
api_key <- "eyJvcmciOiI1YjNjZTM1OTc4NTExMTAwMDFjZjYyNDgiLCJpZCI6ImIzYzBmMDZkMDljNjQ3ZWFhMzM3ZGFlOTVkYTE2MzdkIiwiaCI6Im11cm11cjY0In0="

#Coordenadas para API
coords <- lapply(
  seq_len(nrow(municipios)),
  function(i) c(
    municipios$longitud[i],
    municipios$latitud[i]
  )
)

#Llamar a la Matrix API (Distancia por carretera)

url <- "https://api.openrouteservice.org/v2/matrix/driving-car"

response <- POST(
  url,
  add_headers(
    Authorization = api_key,
    `Content-Type` = "application/json"
  ),
  body = list(
    locations = coords,
    metrics   = list("distance"),
    units     = "km"
  ),
  encode = "json"
)

#Verificar respuesta
response$status_code

#Matriz de distancias
result <- fromJSON(
  content(response, as = "text", encoding = "UTF-8")
)
dist_matrix <- as.matrix(result$distances)
rownames(dist_matrix) <- municipios$Municipio
colnames(dist_matrix) <- municipios$Municipio

#Guardar matriz
write.xlsx(
  dist_matrix,
  file = "distancias_municipios_cordoba.xlsx",
  rowNames = TRUE
)

5. Problema de ruteo de vehículos (VRP: Vehicle Routing Problem).

El VRP es una generalización del TSP. En vez de un solo vehículo, en el VRP se cuentan con múltiples vehículos que salen de un depósito, atienden clientes (vértices) con demanda conocida y regresan al depósito, minimizando el costo total del transporte.

Contextualmente de tiene un depósito central, \(n\) clientes dispersos geográficamente, cada cliente pide cierta cantidad de producto o tiene una demanda y se tienen \(m\) vehículos con capacidad limitada para atender la demanda conocida. Se quiere decidir qué clientes atiende cada vehículo y en qué orden visitarlos minimizando el costo, la distancia, el tiempo.

El modelo base del VRP asume:

1. Un único depósito o nodo de partida.
2. Todos los vehículos son iguales, es decir, la flota es homogénea.
3. Cada cliente se visita una sola vez.
4. Cada vehículo:

• Sale del depósito.
• Visita clientes.
• Regresa al depósito.

5. Existe una capacidad máxima \(Q\).
6. Demanda conocida por cliente \(d_i\).
7. Costos determinísticos.
8. No existen ventanas de tiempo.

• Conjuntos

  • \(V\): conjunto de vértices, \(V: \{0,1,2,..., n \}\). \(0\): depósito, \(1,...,n\): clientes.
  • \(K\): conjunto de vehículos, \(K: \{1,2,...,m \}\).

• Parámetros

  • \(c_{ij}\): costo de recorrer desde \(i \in V\) hasta \(j \in V\).
  • \(d_i\): demanda del cliente \(i \in V, i \ne 0\).
  • \(Q\) capacidad del vehículo.
  • \(m\): número de vehículos, tamaño de flota.

• Variables de decisión

  • Variable binaria de arco o recorrido entre dos vértices.

\[ x_{ij} = \begin{cases} 1, & \text{si algún vehículo } k \in K \text{ viaja de } i \in V \text{ a } j \in V \\ 0, & \text{dlc } \end{cases} \]

• Variable auxiliar para contabilizar la carga (MTZ de capacidad)

\[ u_i: \text{ carga acumulada al llegar al nodo } i \]

Función objetivo

\[ Min~Z: \sum_{i \in V} \sum_{j \in V,~i \ne j} c_{ij}~x_{ij} \]

Restricciones

• Cada cliente se visita una vez

  • Entrada

  \[ \sum_{i \in V,~i \ne j} x_{ij} =1,~~~\forall j \in V~ \backslash ~\{ 0\} \]

  • Salida

  \[ \sum_{j \in V,~i \ne j} x_{ij} =1,~~~\forall i \in V ~ \backslash~ \{ 0\} \]

• Flujo del depósito

  • Número de rutas que salen del depósito

\[ \sum_{j \ne 0} x_{0j} = m \]

• Número de rutas que regresan al depósito

\[ \sum_{i \ne 0} x_{i0} = m \]

• Eliminación de subtours y capacidad

\[ u_i - u_j + Qx_{ij} \leq Q- d_j~\forall i,j \in V ~\backslash~ \{ 0\}, i \ne j \]

• Límites de carga

\[ d_i \leq u_i \leq Q ~~~\forall i \ne 0 \\ \\ u_0 = 0 \]

• Naturaleza de las variables

\[\begin{align} x_{ij} &\in \{0,1 \} \\ u_i &\ge 0 \end{align}\]

Ejemplo 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.

Librerías necesarias

# Cargar librerías
library(ompr)
library(ompr.roi)
library(ROI)
library(ROI.plugin.symphony)
library(readxl)
library(dplyr)
library(leaflet)

Datos

# Dataframe de distancias
distancias_df <- read_excel("distancias_municipios_cordoba.xlsx")

#Vector de municipios
municipios_vec <- distancias_df[[1]]

#Matrix de demandas
dist_matrix <- as.matrix(distancias_df[, -1])

#Nombres de filas y columnas de la matriz
rownames(dist_matrix) <- municipios_vec
colnames(dist_matrix) <- municipios_vec

#Cantidad de nodos
n <- nrow(dist_matrix)

#Datafrma de demandas
demanda_df <- read_excel("demanda.xlsx")

#Orden de demanda
d <- demanda_df[[2]]

#Capacidad por vehículo
Q <- 1.5

#Tamaño de la flota
m <- 12 

# Depósito
nodo_inicial <- which(municipios_vec == "Monteria")

Modelo

modelo <- MIPModel() %>%
  add_variable(x[i, j], i = 1:n, j = 1:n, type = "binary") %>%
  add_variable(u[i], i = 1:n, lb = 0, ub = Q) %>%
  set_objective(sum_expr(dist_matrix[i, j] * x[i, j], i = 1:n, j = 1:n, i != j), "min") %>%
  add_constraint(x[i, i] == 0, i = 1:n) %>%
  add_constraint(sum_expr(x[i, j], i = 1:n) == 1, j = setdiff(1:n, nodo_inicial)) %>%
  add_constraint(sum_expr(x[i, j], j = 1:n) == 1, i = setdiff(1:n, nodo_inicial)) %>%
  add_constraint(sum_expr(x[nodo_inicial, j], j = 1:n) == m) %>%
  add_constraint(sum_expr(x[i, nodo_inicial], i = 1:n) == m) %>%
  add_constraint(u[i] - u[j] + Q * x[i, j] <= Q - d[j], i = setdiff(1:n, nodo_inicial),j = setdiff(1:n, nodo_inicial), i != j) %>%
  add_constraint(u[nodo_inicial] == 0) %>%
  add_constraint(u[i] >= d[i], i = setdiff(1:n, nodo_inicial))

Resultado

resultado <- solve_model(modelo, with_ROI(solver = "symphony",control = list(time_limit = 1200), verbosity = -1))

Extraer arcos

solucion <- get_solution(resultado, x[i, j]) %>% filter(value > 0.5)
print(solucion)
cat("\nDistancia total:", objective_value(resultado), "km\n") #cat imprime texto formateado; \n salto de línea; Distancia total: es el texto; objective_value(resultado) es el valor de la función objetivo; km es la unidad

Reconstruir rutas

routes <- list() #Crea lista vacía para almacenar las rutas finales

salidas <- solucion$j[solucion$i == nodo_inicial] #Busca los nodos destino j tales que existe un arco activo, esto es, x_0j = 1, todos los destinos desde donde sale un vehículo desde Montería.


for(k in seq_along(salidas)){ #Recorre cada salida del depósito
  
  actual <- salidas[k] #El nodo actual es el primer municipio visitado por el vehículo k
  ruta <- c(nodo_inicial, actual) #Inicializa la ruta empezando en el depósito y luego el primer municipio.
  
  while(TRUE){ #Empieza un ciclo infinito que se romperá cuando el vehículo regrese al depósito
    
    siguiente <- solucion %>%
      filter(i == actual) %>%
      pull(j) #Busca un nodo j tal que x_actual,j = 1; encuentra el siguiente municipio visitado desde el nodo actual.
    
    if(length(siguiente) == 0 || siguiente == nodo_inicial) break #Condición de parada, se frena cuando: 1. No existe sucesor, 2. El siguiente nodo es el depósito.
    
    ruta <- c(ruta, siguiente) # Actualiza el nodo actual para continuar recorriendo la cadena.
    actual <- siguiente
  }
  
  ruta <- c(ruta, nodo_inicial) # Cierra la ruta regresando al depósito.
  
  routes[[k]] <- municipios_vec[ruta] # Convierte los índices numéricos en nombres de municipios.
}

routes

Mapa por ruta

#Coordenadas
municipios_df <- read_excel("coordenadas.xlsx")


#Normalizar nombres
municipios_vec <- trimws(toupper(municipios_vec)) #Elimina espacios en blanco y convierto a mayúsculas
municipios_df$Municipio <- trimws(toupper(municipios_df$Municipio))

#Mapa
mapa <- leaflet() %>% addTiles() #Crea objeto interactivo, y agrega mapa base

colores <- rainbow(length(routes)) #Genera una paleta de colores

for(k in seq_along(routes)){ #Itera sobre cada ruta almacenada en routes.
  
  r <- toupper(routes[[k]]) #Extrae la ruta k-ésima y los pone en mayúsculas
  
  ruta_df <- municipios_df %>%
    filter(Municipio %in% r) %>% #Selecciona municipios de la ruta
    slice(match(r, Municipio))   # Reordena las filas según el orden de visita real.
  
  mapa <- mapa %>% #Actualiza el mapa
    addPolylines( #Dibuja lineas de la ruta
      data = ruta_df,
      lng = ~longitud,
      lat = ~latitud,
      color = colores[k], #Usa el color asignado
      weight = 4, #Grosor de la linea
      label = paste("Ruta", k) # Nombre de la ruta al pasar el cursor
    ) %>%
    addCircleMarkers( #Dibuja marcadores circulares para los vértices.
      data = ruta_df,
      lng = ~longitud,
      lat = ~latitud,
      label = ~Municipio,
      radius = 5,
      color = colores[k]
    )
}
mapa