SciELO - Scientific Electronic Library Online

 
vol.84 issue202Probabilistic analysis of the active earth pressure on retaining wall for c -φ soil backfill under seismic loading conditionsProcedure for prioritisation Critical Success Factors author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • On index processCited by Google
  • Have no similar articlesSimilars in SciELO
  • On index processSimilars in Google

Share


DYNA

Print version ISSN 0012-7353

Dyna rev.fac.nac.minas vol.84 no.202 Medellín July/Sept. 2017

https://doi.org/10.15446/dyna.v84n202.62382 

Artículos

Un Algoritmo GRASP híbrido para el 2eCVRP

A hybrid GRASP algorithm for 2eCVRP

Javier Arias-Osorioa 

Andrés Fernando Niño-Sáenza 

1a Escuela de Estudios Industriales y Empresariales, Universidad Industrial de Santander, Bucaramanga, Colombia. jearias@uis.edu.co, andresfernando347@gmail.com


Resumen

Avances recientes en la investigación de problemas de ruteo, han abordado extensiones del clásico Problema de Ruteo de Vehículos (VRP), como lo es el Problema de Ruteo de Vehículos de Dos Escalones (2E-CVRP), en el cual se aborda el diseño de rutas en una cadena de suministro de dos escalones. El primero de ellos que conecta la carga desde un depósito central hasta su consolidación en depósitos intermedios denominados satélites, y el segundo que enlaza la carga de los satélites con el cliente final.

Para la solución del 2E-CVRP se optó por implementar un híbrido metaheurístico, la primera técnica denominada GRASP se enfoca en la formación de una solución inicial; para dar lugar al segundo método, designado como recocido simulado, en el que por medio de los operadores 2-opt, Or-Opt y Exchange, se intensifica la búsqueda de mejora de la solución inicial. Este algoritmo presenta buenos resultados para casos propuestos de la literatura.

Palabras clave: ruteo de vehículos; dos escalones; metaheurísticas; GRASP; recocido simulado.

Abstract

Recent advances in the investigation of routing problems, have allowed to give with extensions of the classic Vehicle Routing Problem (VRP), such as the Two-Echelon Vehicle Routing Problem (2E-CVRP), in which the aim is routes design in two echelon supply chain. The first one that connects the load from a central depot to its consolidation in intermediate deposits called satellites, and the second that links the load of the satellites with the final customer.

For the solution of the 2E-CVRP was opted to implement a metaheuristic hybrid, the first technique called GRASP focuses on the formation of an initial solution; to give rise to the second method designated as simulated annealing, in which the inquisition for improvement of the initial solution is intensified by means of the 2-opt, Or-opt and Exchange operators. This algorithm shows goods results for instances of literature.

Keywords: vehicle routing; two-echelon; metaheuristics; GRASP; simulated annealing.

1. Introducción

La cadena de suministro juega un papel fundamental en la ejecución de la misión de las organizaciones, principalmente a través de la conexión entre proveedores, fabricantes, almacenes y distribuidores, buscando brindar un producto o servicio de calidad al cliente en el momento preciso, al mismo tiempo que se busca optimizar el uso de los recursos disponibles[1].

En este artículo se aborda esa optimización asociada al transporte de mercancías en una cadena de suministro de dos escalones, donde en cada escalón hay un ruteo diferente, pero a la vez complementario.

1.1. Justificación

Las compañías que tienen claro su enfoque de mejora continua, buscan prestar atención a las operaciones logísticas, en pro de ser más eficientes y competitivos en el mercado.

El flujo de los productos se presenta en redes, y específicamente uno de los tipos de redes de mayor relevancia en las cadenas de suministro son las redes de distribución, en las que la demanda de un conjunto de clientes debe ser consolidada en instalaciones fijas para su posterior distribución. La circulación de la demanda de productos fluye por medio de canales de distribución. Estos canales son de vital importancia para el acceso al mercado, y en medio de los fabricantes y el consumidor, generalmente se encuentran los centros de distribución, esto debido a restricciones políticas, físicas o del mercado, entre otras.

Los tipos de redes de distribución por los que puede optar una compañía son diversos y la manera en que se transporta la mercancía puede darse vía terrestre, marítima, aérea o férrea. Sin embargo, un elemento común de las redes de distribución terrestre es el aspecto del ruteo, que considera la manera apropiada de distribuir la mercancía a través de la cantidad de niveles que haya en un determinado sistema.

La investigación de redes de distribución escalonada ha sido abordada recientemente y tiene gran aplicabilidad en “áreas como el servicio de entrega express, distribución de productos comestibles e hipermercados, comercio electrónico y servicios de entrega a casa, distribución de prensa, y logística en la ciudad” [2]; siendo de alto interés el ruteo de este tipo de redes el cual con un apropiado planteamiento, permite mejorar la calidad de vida de las personas, ya que de manera intrínseca considera aspectos como la disminución de la contaminación y el ruido en las ciudades.

Los avances dentro de la ciencia de la computación han permitido descubrir distintos métodos para la solución de problemas complejos como lo son los problemas de ruteo, dentro de los cuales se encuentra el 2E-CVRP (problema de ruteo de vehículos de dos escalones con capacidad limitada); estos métodos de solución apuntan a encontrar buenas soluciones en tiempos oportunos, permitiendo considerar nuevas formas de tratarlos y trabajarlos a nivel teórico.

1.2. Revisión de literatura

El primer desarrollo formal del 2E-CVRP [3] se remonta a 2007, cuando se presenta una formulación del 2E-CVRP basada en el flujo de la mercancía al hacer uso de un algoritmo exacto de ramificación y acotamiento, el cual les permite encontrar soluciones optimas a casos de hasta 20 clientes y dos satélites; además en [2] se introduce la familia de problemas de ruteo de vehículos de dos escalones y se modela la versión básica de esta familia de problemas, es decir el 2E-CVR. Los autores proponen un modelo de programación lineal entera mixta (MILP) en el que las variables de decisión principalmente guardan relación con el flujo de carga en cada arco; además introducen desigualdades válidas para el problema. También, se refieren a los algoritmos de optimización aplicados al problema como heurísticas basadas en matemáticas las cuales combinan elementos de programación matemática con metaheurísticas; estas heurísticas matemáticas se basan en información recuperada de la solución óptima de la relajación lineal del modelo.

En [4] se diseñan dos heurísticas para el 2E-CVRP, que se enfocan en separar el primer escalón del segundo. Ambas tienen un enfoque de dos fases; en la primera fase de ambas heurísticas se encuentra una solución factible para el ruteo del segundo escalón al llevar a cabo una asignación de clientes a los satélites; la diferencia es que una considera el problema del segundo escalón como un conjunto de problemas de ruteo de vehículos con capacidad limitada (CVRP), mientras que la otra lo aborda como un problema de ruteo de vehículo multi-depósito (MDVRP).

En [5] los autores resaltan la aplicabilidad de sistemas de dos escalones en la logística de las ciudades, la cual implica la coordinación de expedidores, transportistas y movimientos, así como la consolidación de la carga que se destinará a los clientes en vehículos amigables con el medio ambiente. El problema que abordan es el 2E-CVRP con dependencia de tiempo, sincronización, multi-depósito, multi-producto, flotas heterogéneas en cada escalón y ventanas de tiempo.

En [6] se presenta una familia de heurísticas multi-inicio, las cuales se basan en separar la transferencia entre: el depósito al satélite, de la entrega del satélite al cliente al resolver iterativamente estos dos subproblemas mediante un algoritmo de búsqueda local basado en el cambio de asignación del cliente al satélite.

En [7] se exponen desigualdades válidas basadas en el problema del agente viajero (TSP) y el CVRP para reforzar el modelo basado en el flujo; además se definen desigualdades válidas basadas en la interacción entre el ruteo y las variables de activación de arcos, así como también otras desigualdades válidas considerando propiedades de conectividad y propias de una solución factible en problemas de ruteo.

En [8] los autores proponen una mejora a la notación de problemas de localización-ruteo, introducida por Laporte [9], en dicha notación la expresión λ/ M 1 /⋅⋅⋅/ 𝑀 λ−1 representa problemas de localización-ruteo multi-escalonados donde λ es el número de etapas o capas, y 𝑀 𝑖 indica el tipo de recorrido entre dos etapas 𝑖 y 𝑖+1. Si entre las etapas 𝑖 y 𝑖+1 la ruta es dedica a un único cliente (viajes hacia y desde un único nodo), entonces 𝑀 𝑖 =𝑅. Si entre las etapas 𝑖 y 𝑖+1 pueden darse rutas multi-cliente (viajes a través de varios vértices), entonces 𝑀 𝑖 =𝑇. Asimismo, sugieren marcar con una línea el indicativo del tipo de recorrido ( 𝑅 o 𝑇 ) si decisiones de localización tienen que ser tomadas en la etapa 𝑖. En este orden de ideas el 2E-CVRP corresponde a un problema3/T/𝑇.

En [10] el autor enfatiza la aplicabilidad del 2E-CVRP en los niveles de planeación estratégica, táctica y operativa; por lo que los métodos para solución de este tipo de problemas deben ser precisos y rápidos, ya que el problema puede necesitar ser solucionado en muchas ocasiones, como parte de un proceso de optimización extenso. Los resultados muestran que el enfoque de dos escalones es aconsejable en ciertas circunstancias frente al de un solo escalón, tal que el 2E-VRP es aplicable a situaciones como entregas de carga en la ciudad.

En [11] se plantea un modelo de optimización simétrico con una formulación MILP para el 2E-CVRP, el cual debido a la simetría provee límites inferiores poco favorables, por lo que los autores proponen una formulación alternativa la cual consiste en una relajación al 2E-CVRP que elimina la simetría y los límites inferiores poco favorables tras lo cual los autores idean una prueba de factibilidad y un esquema de ramificación especializado para obtener soluciones enteras factibles.

En [12] los autores agregan otros componentes del costo, como costos fijos por usar los arcos, costos de operación y costos ambientales; los autores estudian el impacto de la congestión de tráfico en el costo del viaje, por medio de análisis situacionales al variarlos de tal forma que el costo del viaje cambia en diferentes escenarios, pero es constante en cada uno de ellos.

En [13] se expone una heurística de búsqueda local adaptativa y extensa (ALNS) para el 2E-CVRP y el problema de localización y ruteo (LRP), la idea del algoritmo ALNS es llevar a cabo iteraciones en las que se remueve un cliente de la solución actual por medio de un operador destructor, para después, por medio del uso de un operador reparador, reinsertar al cliente en otra posición. El desempeño del algoritmo logra mejorar los resultados en [3].

En [14] se presenta un modelo de programación lineal entera (ILP) para el 2E-CVRP, así como la implementación de dos algoritmos híbridos de ramificación y acotamiento, y generación de columnas; uno de los algoritmos considera rutas que satisfacen la condición primaria, y el otro relaja dicha restricción al fijar precios a las rutas.

En [15] se modela el 2E-CVRP a través de una formulación ILP la cual permite relajaciones enteras y continuas. El algoritmo exacto evaluado descompone el problema en un conjunto limitado de MDVRPs con restricciones laterales.

En [16] se aborda el 2E-CVRP con una heurística que combina el método GRASP (procedimiento de búsqueda voraz adaptativo aleatorio) con un procedimiento de re-encadenamiento de trayectorias (Path Relinking). Al igual que en otros métodos de solución; el 2E-CVRP se divide en dos subproblemas, tras lo cual se efectúan las siguientes fases: primero, un GRASP genera soluciones, segundo, si la solución no es factible se ejecuta una búsqueda de factibilidad, tercero, una búsqueda local que busca mejorar la solución encontrada, y cuarto, un re-encadenamiento de trayectorias.

En [17] se estudia una heurística hibrida de dos fases para el 2E-CVRP la cual combina GRASP con una búsqueda local variable descendente (VND); este algoritmo resuelve el problema de abajo hacia arriba, primero un GRASP genera arreglos de todos los clientes y los asigna a los satélites hasta presentarse una asignación factible, tras lo cual al resolver el problema del primer escalón se obtiene la solución completa a la asignación factible encontrada, posteriormente el método VND explora en las vecindades aplicadas. Acorde con la literatura revisada, esta heurística híbrida genera mejores resultados que las demás.

En [18] se realiza un estudio de la literatura referente a problemas de ruteo de dos escalones. De manera detallada abordan el problema localización-ruteo de dos escalones (2E-LRP), el 2E-VRP, y el problema de ruteo de camiones y remolques (TTRP) al definir en qué consisten, muestran un modelo matemático para cada uno de ellos, describen parte de la literatura más relevante y resumen gran parte de esta información en tablas y gráficos.

En [19] se diseña un algoritmo híbrido el cual acopla un algoritmo voraz, un algoritmo de optimización basado en colonia de hormigas (ACO) y un algoritmo de búsqueda local. A pesar de que el algoritmo genera soluciones de menor calidad comparadas a las soluciones de los algoritmos más precisos, éste tiene un mejor desempeño en cuanto a la velocidad de solución.

En [20] se da a conocer su estudio del 2E-CVRP dependiente del tiempo donde con un modelo MILP, consideran la distancia transitada, la velocidad del vehículo, la carga, las múltiples zonas de tiempo (distintos días en los que los vehículos viajan) y las emisiones. Como conclusión principal, mencionan que un sistema de distribución de dos escalones provee una solución amigable con el medio ambiente, mientras que un sistema de distribución de un único escalón entrega soluciones de menor costo.

En [21] se incorporan restricciones propias de la logística en la ciudad, como lo son restricciones de ventanas de tiempo, restricciones de sincronización y viajes múltiples en el segundo escalón. Para hacer frente al problema, los autores usan una heurística de tipo ALNS la cual proporciona buenas soluciones en tiempos apropiados.

2. Descripción y formulación matemática del problema

El 2E-CVRP consta de dos niveles o escalones, y el modelo matemático, se enfoca en solucionar el ruteo para ambos de forma integral, a la vez que se considera minimizar el costo del problema; se entiende como ruta de primer nivel a una que conecta un depósito central con uno o más centros de distribución conocidos como satélites, y como ruta de segundo nivel a una que conecta a un satélite con uno o más clientes que requieren una demanda, las rutas del segundo escalón requieren ser asignadas a alguno de los satélites, mientras que las del primer escalón no. Se cuenta con dos flotas homogéneas de vehículos, una para cada nivel, siendo la flota del segundo nivel menor en capacidad respecto a la del primer nivel. Usualmente, restricciones de capacidad son consideradas para los satélites y vehículos. La demanda de cada cliente no puede ser dividida en los vehículos del segundo nivel; sin embargo, los vehículos del primer nivel pueden transportar la carga de uno o más clientes, así como servir más de un satélite en la misma ruta; cabe aclarar que los satélites pueden ser visitados por distintos vehículos y las entregas de estos pueden ser divididas. Finalmente, la carga debe ser consolidada en los satélites para su posterior envío a los clientes. En la Fig. 1 se evidencia una representación gráfica del 2E-CVRP:

Fuente: Adaptada de [17].

Figura 1 Ilustración del 2E-CVRP. 

De acuerdo con las especificaciones planteadas en [17], dado un grafo ponderado, completo y no dirigido 𝐺=(𝑉,𝐴) en el que el conjunto de vértices 𝑉 se conforma por: un depósito 𝑉 0 , un subconjunto 𝑉 𝑆 y un subconjunto 𝑉 𝐶 ; dichos subconjuntos con 𝑛 𝑠 satélites y 𝑛 𝑐 clientes, respectivamente.

Respecto al conjunto de aristas 𝐴, se compone por los arcos (𝑖,𝑗) que representan el camino que conecta a los vértices al cual se le asocia un costo por su recorrido 𝐶 𝑖𝑗 ; en el que estos costos satisfacen la desigualdad triangular en los dos escalones. A cada uno de los clientes 𝑖∈𝑉 𝐶 se le asocia una demanda 𝑑 𝑖 . Los vehículos tienen una capacidad limitada 𝑄 1 y 𝑄 2 , para el nivel 1 y 2 respectivamente. El número total de vehículos en el primer escalón será de 𝑚 1 y para el segundo de 𝑚 2 . Además, cada satélite 𝑘∈𝑉 𝑠 tiene una capacidad 𝐵 𝑘 la cual restringe la cantidad de demanda que llega a este desde las rutas de primer nivel. Como última limitante se tiene que 𝑚 𝑘 vehículos del primer nivel están disponibles en cada satélite 𝑘∈𝑉 𝑠 .

La función objetivo del problema se expresa como la minimización de:

En la que se busca minimizar los costos de recorrido asociados a las rutas de primer y segundo nivel.

Respecto a las restricciones se tiene que:

  • Cada cliente debe ser visitado por un solo vehículo el cual efectúa una única ruta de segundo nivel.

  • En cada satélite, el número de rutas posibles no puede ser mayor al de la cantidad de vehículos disponibles para efectuarlas.

  • El número total de rutas activas en el segundo escalón no puede exceder al número total de vehículos en este.

  • La capacidad de cada satélite debe ser mayor o igual a la demanda que se consolidará en este, para su respectivo envío a los clientes.

  • El número total de rutas activas en el primer escalón no puede exceder al número total de vehículos en este.

  • La cantidad de demanda que llega desde las rutas de primer nivel debe ser igual a la cantidad demandada por los clientes asignados en cada satélite.

  • Para cada una de las rutas posibles en las que participa un vehículo de primer nivel, la capacidad de este no puede ser excedida por la demanda de los satélites a los que sirve.

  • Las variables del modelo son de tipo binario y entero.

3. Metaheurística híbrida aplicada al 2E-CVRP

En la primera fase del híbrido con el que se abordó el 2E-CVRP se aplicó un procedimiento de búsqueda voraz adaptativo aleatorio (GRASP), el cual se adaptó del expuesto en [18], dentro de esta fase se emplea un procedimiento de división (spliting procedure) basado en Prins [22]; este procedimiento es vital para la construcción de una solución inicial en el marco de esta fase.

La segunda fase comprende los esfuerzos por mejorar las soluciones encontradas en la primera fase, para ello se utiliza la metaheurística de recocido simulado (SA) dentro de la cual se establecen las vecindades de búsqueda pertinentes.

3.1. Procedimiento de búsqueda voraz adaptivo

Este método iterativo consiste en generar de forma aleatoria recorridos [17] como los del problema del agente viajero (TSP), estos recorridos que forman parte del segundo escalón, cubren la cantidad total de clientes y son la base para adaptarles el procedimiento de división, el cual permite romper el recorrido TSP en viajes asignados a cada uno de los satélites. Al conocer los clientes de cada viaje y su asignación a su respectivo satélite, es posible encontrar la demanda de cada satélite, la cual equivaldría a la suma de las demandas de los clientes asignados al mismo.

Con la demanda de los satélites se procede para el primer escalón a hacer tantos envíos de carga completa como sea posible desde el depósito a cada uno de los satélites, cuando la demanda remanente de cada satélite es menor a la capacidad del vehículo de primer escalón se da la condición necesaria para aplicar nuevamente el proceso de división, el cual permitirá dar con los viajes que servirán la demanda faltante de los satélites.

3.1.1. Procedimiento de división

Para una determinada permutación 𝑇 de clientes, se procede a crear viajes a los clientes, asumiendo que un vehículo contara con capacidad infinita y pudiera visitar a todos los clientes de acuerdo al orden establecido en esta, así como a realizar asignaciones de estos viajes a los satélites.

Este procedimiento [23] funciona por medio de un grafo auxiliar 𝐻=(𝑋,𝐴,𝑍) en el que:

  • X representa el conjunto de nodos desde 0 hasta el número de clientes.

  • A simboliza el conjunto de los arcos 𝑖,𝑗 , 𝑖<𝑗, que representan los viajes factibles, es decir los viajes que cumplen con la siguiente restricción de carga:

  • La suma de las demandas de los clientes 𝑝 𝑖+1 a 𝑝 𝑗 visitados en el viaje no exceda la capacidad del vehículo:

  • Z representa el conjunto de los costos de los viajes 𝑧 𝑖𝑗 de los arcos 𝑖,𝑗 .

A modo de ejemplo [24], dada una permutación =[ 𝑝 3 , 𝑝 5 , 𝑝 1 𝑝 2 , 𝑝 4 ] , y siendo este un vector que contiene el arreglo aleatorio de cinco clientes, se requiere además conocer la capacidad del vehículo de segundo escalón, la cual será 𝑄 2 =15, las demandas de los clientes se encuentran entre paréntesis en la Fig. 2, la cual además revela los costos (asumidos como distancia) de ir de un nodo a otro.

Fuente: Elaboración propia.

Figura 2 Grafo de la permutación T. 

Nótese la posibilidad de los futuros cortes (viajes) de ser asignados a los satélites 𝑆 2 o 𝑆 1 . El detalle del grafo auxiliar 𝐻 del ejemplo se evidencia en la Fig. 3:

Fuente: Elaboración propia.

Figura 3 Grafo auxiliar del procedimiento de división. 

Cada uno de los arcos de la Fig. 3, representa un viaje asignado al satélite de menor costo, para el caso del arco 𝑝 1 𝑝 2 𝑝 4 el valor del viaje desde el satélite uno seria de 125 (resultado de 35+30+20+40), mientras que desde el satélite dos el valor sería de 135 (resultado de 50+30+20+35); por lo cual la asignación de menor costo corresponde al satélite uno. De tal manera que el camino más corto hasta el último nodo del grafo auxiliar, el cual toma un valor de 195, indica cual es la mejor manera de dividir la permutación, que para éste caso se encuentra demarcado por los dos arcos con línea acentuada, además se evidencia como si se hubiera dividido de otra forma, los costos hubieran aumentado como en el caso de los caminos que dan un valor de 225 y 245. Por lo que el resultado gráfico del procedimiento de división aplicado al ejemplo se evidencia en la Fig. 4.

Fuente: Elaboración propia.

Figura 4 Solución por medio del procedimiento de división. 

3.2. Recocido simulado

La técnica recocido simulado [25], se basa en un concepto propio de la ciencia de materiales, denominado recocido total [24] y consiste en una solución actual 𝑠, en la que el costo asociado a esta se compara con una solución 𝑡∈𝑁(𝑠), por lo cual la solución 𝑡 pertenece a la vecindad de la solución actual 𝑠, siendo 𝑠 para el caso de la primera iteración, equivalente a la primera solución obtenida de la metaheurística GRASP; de tal manera que si el costo de la solución 𝑡 es menor o igual al costo de la solución actual 𝑠, entonces 𝑡 tomaría el lugar de 𝑠; para el caso contrario en el que el costo de la solución 𝑡 sea mayor al de 𝑠, se le asigna una probabilidad a 𝑡 de llegar a ser la solución actual, dicha probabilidad junto con la expresión matemática del caso se pueden expresar de la siguiente forma [26]:

Donde 𝑓 𝑠 y ?? 𝑡 son los valores del costo asociados a las soluciones 𝑠 y 𝑡 respectivamente; y 𝑡𝑒𝑚 𝑘 representa la temperatura del sistema en la iteración 𝑘. Además, al analizar la expresión de probabilidad se reconoce como entre más alta sea la temperatura, mayor será la probabilidad de aceptar una solución que no mejora la actual, como la sustituta de 𝑠.

La forma clásica de programar el enfriamiento se presenta en la siguiente expresión:

Por lo cual la velocidad de enfriamiento depende del parámetro 𝛽, en el que valores altos de este, significa que el sistema se enfría de una manera lenta, comparado con valores bajos de 𝛽, los cuales implican un enfriamiento más rápido.

Dado que GRASP garantiza la exploración adecuada del problema, al incrustar el recocido simulado como enfoque de la fase de intensificación, se requiere limitar la exploración de esta, situación que se logra al optar por una temperatura inicial tibia [26].

En el marco de aplicación del recocido simulado, se hace vital la construcción de las vecindades [27], sobre las cuales se llevará a cabo la intensificación de la búsqueda en pro de mejorar la solución de entrada.

3.2.1. Operador de construcción 2-opt

Para una determinada ruta de segundo nivel, al usar 2-opt se remueven dos aristas no adyacentes de la ruta, y se posicionan dos nuevas de tal forma que la dirección en la cual se atendían los clientes dentro de los dos arcos, se da al revés [29].

A modo de ejemplo, en la Fig. 5, se evidencia como los arcos que conectan a los clientes ( 𝑝 11 , 𝑝 7 ) y ( 𝑝 4 , 𝑝 9 ) son expelidos; para dar lugar a los arcos ( 𝑝 11 , 𝑝 4 ) y ( 𝑝 7 , 𝑝 9 ), los cuales hacen parte de la nueva ruta, en la cual además, se ve claramente como el sentido de los arcos ( 𝑝 4 , 𝑝 2 ), ( 𝑝 2 , 𝑝 1 ), ( 𝑝 1 , 𝑝 5 ), ( 𝑝 5 , 𝑝 8 ) y ( 𝑝 8 , 𝑝 7 ) se dan en dirección contraria respecto a los arcos de la ruta no modificada.

Fuente: Elaboración propia.

Figura 5 Aplicación del operador 2-opt. 

3.2.2. Operador de construcción Or-opt

Este operador [30] reubica la posición de un nodo cliente perteneciente a una ruta de segundo nivel, de tal forma que en el gráfico de la Fig. 6, se expone como la eliminación de las aristas ( 𝑝 5 , 𝑝 2 ), ( 𝑝 2 , 𝑝 7 ) y ( 𝑝 3 , 𝑝 6 ), las cuales forman parte de la ruta inicial; y la disposición de tres nuevos arcos, es decir, las aristas ( 𝑝 5 , 𝑝 7 ), ( 𝑝 3 , 𝑝 2 ) y ( 𝑝 2 , 𝑝 6 ), dan como resultado un cambio de posición del cliente 𝑝 2 en la ruta final modificada.

Fuente: Elaboración propia.

Figura 6 Aplicación del operador Or-opt. 

3.2.3. Operador de construcción Exchange

En este caso [30] se troca la ubicación de dos nodos clientes de una ruta de segundo nivel, los cuales no tienen limitante frente a su adyacencia; por lo cual en el primer grafo de la Fig. 7, al suprimir los arcos (p_5,p_2), (p_2,p_7), (p_3,p_6) y (p_6,p_1), y dar cabida a los arcos (p_5,p_6), (p_6,p_7), (p_3,p_2) y (p_2,p_1), se suscita el intercambio de los clientes p_2 y p_6.

Fuente: Elaboración propia.

Figura 7 Aplicación del operador Exchange. 

3. Extensión del método híbrido

La forma en la que funciona el procedimiento de división lleva a inferir que una permutación construida con base en criterios que consideren la distancia entre los clientes tiene una posibilidad alta de dar con soluciones iniciales de un costo bajo.

La extensión del método híbrido consiste en generar una matriz de tamaño 𝑛∗𝑛, con 𝑛 igual al número de clientes del problema; cada fila de la matriz corresponde a una permutación, y cada una de estas permutaciones debe empezar por un cliente 𝑖 una única vez. A razón de que se cuenta con las coordenadas de los clientes, se procede a comparar la distancia euclidiana entre el cliente más reciente asignado a la permutación y todos los demás clientes que no forman parte de esta, de tal manera que el cliente con la menor distancia se asigna en la posición 𝑗+1.

Una vez se cuenta con la matriz de permutaciones, a cada una de ellas se le aplica el procedimiento de división; en el primer escalón se crea una segunda matriz de permutaciones que a diferencia de la matriz del segundo escalón, contiene todas las permutaciones posibles del conjunto de satélites, la razón por la que se generan todas las permutaciones en la matriz del primer escalón, es su factibilidad computacional en el momento de solucionar el problema, debido a la menor cantidad de satélites respecto de clientes, después se resuelve el problema para tantas permutaciones como clientes hay en el problema (matriz del segundo escalón) y cada una de estas permutaciones evaluadas con todos los ordenamientos posibles de los satélites que se solucionan por medio del procedimiento de división en el primer escalón.

A modo de ejemplo, para la construcción de la matriz de segundo escalón se considera el conjunto de clientes 𝐶=[ 𝑝 1 , 𝑝 2 , 𝑝 3 𝑝 4 , 𝑝 5 , 𝑝 6 ], para el que se debe generar seis permutaciones, en la que en cada una se inicie con un cliente diferente; como se ve en la Fig. 8, la primera permutación comienza con 𝑝 1 , el cual se ubica en la posición (1,1) de la matriz (ver Fig. 11).

Fuente: Elaboración propia.

Figura 8 Asignación inicial de la primera permutación.  

Con los valores correspondientes a las coordenadas de los clientes de la Tabla 1, se procede a calcular la distancia euclidiana entre 𝑝 1 y los demás clientes que aún no se asignan en la permutación, es decir ( 𝑝 2 , 𝑝 3 𝑝 4 , 𝑝 5 , 𝑝 6 ).

Tabla 1 Coordenadas de los clientes. 

Fuente: Elaboración propia.

En la Tabla 2, se presentan los valores correspondientes a la distancia entre 𝑝 1 (el cliente inicial de la primera permutación) y los clientes ( 𝑝 2 , 𝑝 3 𝑝 4 , 𝑝 5 , 𝑝 6 ).

Tabla 2 Comparación de distancias. 

Fuente: Elaboración propia.

Como se refleja en la Tabla 2, el cliente más cercano a 𝑝 1 es 𝑝 2 , puesto que la distancia 𝑝 1 , 𝑝 2 es igual a 8,544, siendo este el escalar de menor cuantía entre las demás distancias, debido a ello, 𝑝 2 toma la posición (1,2) dentro de la matriz (ver Fig. 11), y por tanto la segunda posición de la primera permutación como se detalla en la Fig. 9.

Fuente: Elaboración propia.

Figura 9 Segunda asignación de la primera permutación. 

Este procedimiento se aplica de igual forma para 𝑝 2 , se resalta que en esta segunda iteración, no se debe considerar 𝑝 1 en la comparación de distancias con 𝑝 2 , puesto que ya se ubicó dentro de la permutación en construcción. Al terminar las repeticiones de la primera permutación se obtiene la primera fila de la matriz la cual se muestra en la Fig. 10.

Fuente: Elaboración propia.

Figura 10 Primera permutación de la matriz. 

Ya completa la primera fila de la matriz del segundo escalón, se replica el procedimiento aplicado hasta completarla, teniendo presente que el primer cliente de cada fila no debe ser igual en ninguna de estas, en la Fig. 11, se muestra la matriz final del ejemplo.

Fuente: Elaboración propia.

Figura 11 Matriz de permutaciones completa. 

Finalmente, se considera crear la segunda matriz de permutaciones respecto de un determinado conjunto de satélites, sobre la cual se prueba cada una de las permutaciones de la primera matriz.

4. Resultados de la experimentación

Para verificar los algoritmos desarrollados, se tomaron como referencia los casos propuestos en [2]. De estas se seleccionaron diez casos, las cuales se probaron por medio de los algoritmos de las técnicas desarrollados en MATLAB®.

En [13] los autores elaboran un resumen de los conjuntos de casos pertinentes al 2E-CVRP el cual se presenta en la Tabla 3, donde se denotan la primera columna S(Conjuntos) y la segunda I (casos). En el caso de la primera fila, seis casos que forman parte del conjunto número dos, tienen como parámetros los valores propios de la fila en la que están presentes. El número de satélites se indican en la columna 𝑚, la columna 𝑛 muestra el número de clientes, las columnas 𝑚 1 y 𝑚 2 presentan el número de vehículos disponibles para el ruteo del primer y segundo escalón respectivamente, y las capacidades de los vehículos del primer y segundo escalón vienen dadas por 𝐾 1 y 𝐾 2 , respectivamente.

Tabla 3 Características de los casos. 

Fuente: Adaptado de [14].

En la Tabla 4, se presentan los resultados correspondientes al valor de la distancia (costo) de las soluciones correspondientes a los casos especificados en la columna I, por medio del método híbrido desarrollado (GRASP+SA) y la primera etapa de éste mismo método, es decir, únicamente haciendo uso de GRASP.

Tabla 4 Resultados GRASP vs. GRASP+SA. 

Fuente: Elaboración propia.

Respecto a la Tabla 4:

La columna nombre indica la manera en la que se le llama al caso en la literatura.

Las columnas GRASP y GRASP+SA contienen los resultados de la aplicación de la primera parte del método híbrido y de la totalidad del mismo respectivamente.

En la última columna se evidencia el porcentaje de diferencia.

De la Fig. 12, se identifica como la etapa de intensificación con respecto de las soluciones iniciales arrojadas por GRASP tiene un efecto positivo en los resultados para todos los casos, de tal manera que las vecindades de búsqueda construidas a partir de los operadores 2-opt, Or-opt y Exchange, dan con soluciones mejor organizadas en cuanto al ruteo. El efecto de la etapa de intensificación se evidencia más profundamente para el caso E-n51-k5-s2-4-17-46, en la que se llega a una variación del -32,01%.

Fuente: Elaboración propia.

Figura 12 Gráfico comparativo GRASP vs. GRASP+SA. 

Respecto a los resultados de la Tabla 5, se considera la comparación entre GRAS+SA y la ampliación de dicha técnica (ver numeral 3.3), en la que varía la forma como se generan las permutaciones iniciales previas al procedimiento de división, dejando a un lado el método aleatorio, lo cual trae efectos positivos sobre las soluciones obtenidas de dicha forma.

Tabla 5 Resultados GRASP+SA vs. GRASP+SA extendido. 

Fuente: Elaboración propia.

Respecto a la Tabla 5, las columnas abarcan el mismo significado que el de la Tabla 4, a excepción de la columna GRASP+SA extendido la cual contienen los resultados de la aplicación de la extensión del método híbrido desarrollado en el que se considera la técnica de construcción de las permutaciones del problema, respectivamente.

Tal como se evidencia en la Tabla 5, la extensión del híbrido permite llegar a resultados de hasta un -48,47% de mejora respecto del GRASP+SA como lo es para el caso “E-n51-k5-s13-19”, lo cual demuestra como la formación de las permutaciones acorde a la extensión favorece la aplicación del procedimiento de división, puesto que la dirección en la que este último construye los viajes o cortes se basa en la distancia, y una permutación en la que se procuren vecinos cercanos dentro de ella permite dar con viajes más cortos. Un panorama general del desempeño de la extensión se muestra en la Fig. 13, donde el rendimiento de esta sobrepasa a todos los casos evaluados, con las diferencias porcentuales más altas en los casos que cuentan con los valores elevados en cuanto a la cantidad de clientes y satélites, esto se debe a que al aumento del dominio de soluciones del problema dificulta la consecución de permutaciones aleatorias de calidad en el GRASP+SA.

Fuente: Elaboración propia.

Figura 13 Gráfico comparativo GRASP+SA vs. GRASP+SA extendido. 

Por último, se efectúa la comparación entre los costos de las soluciones del método GRASP+VND [17] y la técnica propuesta en el presente trabajo, es decir, el GRASP+SA extendido.

Respecto de la Tabla 6, las columnas abarcan el mismo significado que el de la Tabla 5, a excepción de la columna GRASP+VND donde se detallan los valores de las soluciones de [18].

Tabla 6 Resultados GRASP+VND vs. GRASP+SA extendido. 

Fuente: Elaboración propia.

Visualmente la Fig. 14 evidencia similitud entre la extensión del GRASP+SA y los resultados del GRASP+VND, de hecho, en la Tabla 6 se verifica como en el caso “E-n51-k5-s13-19” se reduce en un 2,06% el mejor valor encontrado hasta el momento.

Fuente: Elaboración propia.

Figura 14 Gráfico comparativo GRASP+VND vs. GRASP+SA extendido. 

5. Conclusiones

Acorde con la literatura, se evidencia una tendencia de adicionar nuevas características que modelan situaciones de interés, como las relacionadas con las emisiones al ambiente, ventanas de tiempo, flotas heterogéneas en cada escalón y el impacto de la congestión en los costos, entre otros. De tal forma que se percibe una propensión por modelar nuevos elementos como demandas estocásticas, clientes estocásticos, horizontes de planeación de las rutas y envíos al depósito a través de los satélites.

De acuerdo a la investigación, se evidencia como el 2E-CVRP puede ser configurado tanto a través de un modelo basado en el flujo, así como por medio de un modelo sustentado en la partición de conjuntos, siendo este último más compacto que el primero, manteniendo todas las características propias del problema.

Respecto del método para la solución del 2E-CVRP, se observa el procedimiento de división que toma protagonismo por su función en la construcción de soluciones iniciales en la primera etapa, mientras que, en la segunda son los operadores de construcción de vecindades los que determinan que se intensifique la búsqueda de soluciones de mayor calidad.

En las pruebas del híbrido inicial se refleja cómo la intensificación surte un efecto relevante, puesto que en la primera fase hay un alto componente aleatorio en la generación de las permutaciones dentro de un dominio de posibilidades extenso, sin embargo, con la presentación de la extensión, se percibe que en la etapa de intensificación no impacta sustancialmente la función objetivo, debido a la calidad de las soluciones iniciales.

A partir de lo anterior, se pudo evidenciar en la aplicación de las técnicas metaheurísticas que, el híbrido GRASP+SA funcionó mejor en el 100% de los casos que la técnica GRASP (Tabla 4). También que al incorporar al híbrido las heurísticas de intercambio, en el 100% de los casos dió mejores resultados que el anterior. Y ya comparando el GRASP+SA extendido (método abordado en este artículo) contra la técnica GRASP+VND revisada de la literatura, se muestra (en la Tabla 6) que en 2 de los casos se obtienen resultados con una diferencia por encima del 20%, que en 7 de los casos se obtienen resultados con una diferencia inferior al 15%. Y en uno de los casos se logra mejorar la solución establecida en la literatura con la técnica utilizada, obteniendo una solución 2,06% mejor.

6. Recomendaciones

En el presente trabajo se emplean operadores para la formación de vecindades, estos realizan variaciones internas en las rutas; además de este tipo de operadores, otros se encargan de realizar movimientos entre rutas como se puede evidenciar en la literatura, por lo cual se sugiere analizar los efectos positivos de una u otra clase de operadores, de los cuales se evidencia como el considerar ambos tipos de operadores permite una mayor intensificación en esta etapa.

En el estudio de sistemas de N-escalones [31] se percibe la necesidad de ajustar el modelo del problema localización-ruteo de n-escalones al problema de ruteo de vehículos de n-escalones con capacidad limitada [31], además de crear casos para problemas de más de dos escalones.

Referencias

[1] Fahimnia, B., Luong, L. and Marian, R. An integrated model for the optimization of a two-echelon supply network. Journal of Achievements in Materials and Manufacturing Engineering, 31(2), pp. 477-484, 2008. [ Links ]

[2] Perboli, G., Tadei, R. and Vigo, D., The two-echelon capacitated vehicle routing problem: Models and math-based heuristics. Transportation Science, 45(3), pp. 364-380, 2011. DOI: 10.1287/trsc.1110.0368 [ Links ]

[3] Gonzales, J., Perboli, G., Tadei R. and Vigo, D., The two-echelon capacitated vehicle routing problem. Technical report DEIS OR.INGCE2007/2(R). [ Links ]

[4] Crainic, T., Mancini, S., Perboli, G. and Tadei, R., Clustering-based heuristics for the two-echelon vehicle routing problem. Technical Report Cirrelt-2008-46. [ Links ]

[5] Crainic, T., Ricciardi, N. and Storchi, G., Models for evaluating and planning city logistics systems. Transportation Science , 43, pp. 432-454, 2009. DOI: 10.1287/trsc.1090.0279 [ Links ]

[6] Crainic, T., Mancini, S., Perboli, G. and Tadei, R., Multi-start heuristics for the two-echelon vehicle routing problem. Evolutionary computation in combinatorial optimization: Lecture notes in Computer Science, 6622, pp.179-190, 2011. DOI: 10.1007/978-3-642-20364-0_16 [ Links ]

[7] Perboli, G. and Tadei, R., New families of valid inequalities for the two-echelon vehicle routing problem. Electronic notes in discrete mathematics, 36(c), pp. 639-646, 2010. DOI: 10.1016/j.endm.2010.05.081 [ Links ]

[8] Crainic, T., Sforza, A. and Sterle, C., Location-routing models for two-echelon freight distribution system design. Technical Report Cirrelt-2011-40. [ Links ]

[9] Laporte, G., Location-routing problems, in: Golden, B. and Assad, A., Vehicle routing: Methods and studies, Elsevier Science Publishers, North-Holland: Amsterdam, 1988. pp. 163-198. [ Links ]

[10] Mancini, S., The two-echelon vehicle routing problem. 4or: A Quarterly Journal of Operations Research, 10(4), pp. 391-392, 2012. [ Links ]

[11] Jepsen, M., Spoorendonk, S. and Ropke, S., A branch-and-cut algorithm for The two-echelon capacitated vehicle routing problem. Transportation Science , 47(1), pp. 23-37, 2012. DOI: 10.1287/trsc.1110.0399 [ Links ]

[12] Crainic, T., Mancini, S., Perboli, G. and Tadei, R., Impact of generalized travel costs on satellite location in the two-echelon vehicle routing problem. Procedia - Social and Behavioral Sciences, 39, pp. 195-204, 2012. DOI: 10.1016/j.sbspro.2012.03.101 [ Links ]

[13] Hemmelmayr, V., Cordeau, J. and Crainic, T., An adaptive large neighborhood search heuristic for two-echelon vehicle routing problems arising in city logistics. Computers & Operations Research, 39(12), pp. 3215-3228, 2012. DOI: 10.1016/j.cor.2012.04.007 [ Links ]

[14] Santos, F., Da cunha, A. and Mateus, G., Branch-and-price algorithms for The two-echelon capacitated vehicle routing problem. Optimization Letters, 7(7), pp. 1537-1547, 2012. DOI: 10.1007/s11590-012-0568-3 [ Links ]

[15] Baldacci, R., Mingozzi A., Roberti, R. and Calvo, R. , An exact algorithm for The two-echelon capacitated vehicle routing problem. Operations Research, 61(2), pp. 298-314, 2013. DOI: 10.1287/opre.1120.1153 [ Links ]

[16] Crainic, T., Mancini, S., Perboli, G. and Tadei, R., Grasp with path relinking for the two-echelon vehicle routing problem. Advances in metaheuristics: Operations research/computer science interfaces series, 53, pp. 113-125, 2013. DOI: 10.1007/978-1-4614-6322-1_7 [ Links ]

[17] Zeng, Z., Xu, W. and Xu, Z., A two-phase hybrid heuristic for the two-echelon vehicle routing problem. Chinese Automation Congress (CAC), pp. 625-630, 2013. DOI: 10.1109/CAC.2013.6775811 [ Links ]

[18] Cuda, R., Guastaroba, G. and Speranza, M., A survey on two-echelon routing problems. Computers and Operations Research , 55, pp. 185-199, 2014. DOI: 10.1016/j.cor.2014.06.008 [ Links ]

[19] Jiang, H., An hybrid heuristic algorithm for the two-echelon vehicle routing problem. Information Science and Control Engineering, 1, pp. 1-5, 2014. [ Links ]

[20] Soysal, M., Bloemhof, J. and Bektas, T., The time-dependent two-echelon capacitated vehicle routing problem with environmental considerations. International Journal of Production Economics, 164, pp. 366-378, 2014. DOI: 10.1016/j.ijpe.2014.11.016 [ Links ]

[21] Grangier, P., Gendreau, M., Lehuédé F. and Rousseau, L., An adaptive large neighborhood search for the two-echelon multiple-trip vehicle routing problem with satellite synchronization. European Journal of Operational Research, 254(1), pp. 80-91, 2016. DOI: 10.1016/j.ejor.2016.03.040 [ Links ]

[22] Prins, C., A simple and effective evolutionary algorithm for the vehicle routing problem. Computers & Operations Research, 31, pp. 1985-2002, 2004. DOI: 10.1016/S0305-0548(03)00158-8 [ Links ]

[23] Nguyen, V., Prins, C. and Prodhon, C., Solving the two-echelon location routing problem by a grasp reinforced by a learning process and path relinking. European Journal of Operational Research , 216(1), pp. 113-126, 2012. DOI: 10.1016/j.ejor.2011.07.030 [ Links ]

[24] Avner, S., Introducción a la metalurgia física. 2° ed. México: Mc graw Hill. 1995. [ Links ]

[25] Festa, P. and Resende, M., Hybrid grasp heuristics, in: Abraham, A., et al. Foundations of Computational Intelligence, vol. 3, 2009. pp. 75-100. DOI: 10.1007/978-3-642-01085-9_4 [ Links ]

[26] Glover, F. and Kochenberger, G., Handbook of metaheuristics. 1° ed. New York: Kluwer academic publishers. 2003. [ Links ]

[27] Ahuja, R., Ergun, O., Orlin, J. and Punnen, A., A survey of very large-scale neighborhood search techniques. Discrete Applied Mathematics, 123, pp. 75-102, 2002. DOI: 10.1016/S0166-218X(01)00338-9 [ Links ]

[28] Chen, P., Huang, H. and Dong, X., Iterated variable neighborhood descent algorithm for the capacitated vehicle routing problem. Expert Systems with Applications, 37, pp. 1620-1627, 2010. DOI: 10.1016/j.eswa.2009.06.047 [ Links ]

[29] Potvin, J. and Rousseau, J., An exchange heuristic for routeing problems with time windows. Journal of the Operational Research Society, 46, pp. 1433-1446, 1995. DOI: 10.1057/jors.1995.204 [ Links ]

[30] Braysy, O. and Gendreau, M., Vehicle routing problem with time windows, Part I: Route construction and local search algorithms. Transportation Science , 39(1), pp. 104-118, 2005. DOI: 10.1287/trsc.1030.0057, 10.1287/trsc.1030.0056 [ Links ]

[31] Gonzales, J., The N-echelon location routing problem: Concepts and methods for tactical and operational planning. Technical report halshs-00422492v2. [ Links ]

1How to cite: Arias-Osorio, J.A. and Niño-Sáenz, A.F., Un algoritmo GRASP híbrido para el 2eCVRP. DYNA, 84(202), pp. 16-25, September, 2017.

Recibido: 01 de Febrero de 2017; Revisado: 08 de Junio de 2017; Aprobado: 14 de Junio de 2017

Creative Commons License Este es un artículo publicado en acceso abierto bajo una licencia Creative Commons