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.

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:

Figura 1 Fuente: Adaptada de [17].

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.

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.

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

Figura 8 Fuente: Elaboración propia.

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

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

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.

Figura 9 Fuente: Elaboración propia.

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.

Figura 10 Fuente: Elaboración propia.

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.

Figura 11 Fuente: Elaboración propia.

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.

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

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

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

Figura 12 Fuente: Elaboración propia.

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

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.

Figura 13 Fuente: Elaboración propia.

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

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.

Figura 14 Fuente: Elaboración propia.

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.

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.