1. Introducción

En los últimos años, la complejidad en la distribución de la última milla ha incrementado significativamente, debido al aumento en las ventas de productos que se entregan puerta a puerta; se estima que los envíos generalmente oscilan de 50 a 150 paradas por día [1]. Adicionalmente, según datos del Estudio del E-Commerce 2020 tuvo un balance positivo. Por ejemplo, según una encuesta de Ernst & Young, 63% de las tiendas online calcularon que su crecimiento sería mayor al 10% lo que resulta fundamental elegir un comercio electrónico y que el plazo de entrega propuesto sea rápido. Así mismo, entre los aspectos que denotan un crecimiento, destaca la intención de internacionalización, según el observatorio e-commerce, el 51% de los comercios encuestados prevén internacionalizar su negocio [2]. Dicha actividad requiere un conjunto de acciones propias o ejecutadas por la distribución de última milla, operación que implica el movimiento de la mercancía o producto desde el último punto de distribución hasta el punto de entrega del usuario o destinatario final. En este sentido, se comprende que los inconvenientes no solo tienen que ver con el tráfico, dificultad en zonas peatonales, la reglamentación y regulación de dichas descargas, sino que también se despachan productos de tamaño pequeño a clientes individuales; por ello, este proceso tiene un costo considerable y una eficiencia baja. Además, desencadena disputas corporativas, causando temas relacionados con la contaminación ambiental, la conciencia ecológica, circunstancias laborales en los repartidores, etcétera. Por lo anterior, las empresas y los investigadores han centrado su atención en determinar rutas con el mínimo costo posible que permitan satisfacer a la totalidad de sus clientes y asignar los mejores recursos para maximizar los niveles de servicio y minimizar los costos en el segmento de transporte final [3]. Encontrado, por ejemplo, que los costos de transacción y transporte juntos suman entre 20% (para los agricultores menos remotos) y más del 50% (para los agricultores más remotos) a los precios de los fertilizantes cobrados en el centro de distribución de insumos [4].

El uso de Vehículos Aéreos No Tripulados o más conocidos como drones, ha ido incrementando considerablemente debido a los avances tecnológicos que han permitido mejorar sus características técnicas además de tener ventajas sobre los medios de transporte terrestres, una de ellas es el bajo impacto negativo que tiene sobre el medio ambiente, la eliminación de atascos de tráfico, entre otras; todo esto ha motivado a usarlos en una gran variedad de aplicaciones. Entre las primeras aplicaciones se tiene las operaciones militares, donde se han usado en labores de vigilancia, desminado de territorios y bombardeos; otra aplicación ha sido en operaciones de logística humanitaria, entre las que se destacan el mapeo de zonas de desastre, apoyo en operaciones de rescate y la distribución de insumos médicos [5]. Viendo estos avances y las ventajas que podría generar, el CEO de Amazon, anunció en 2013 el Prime air, un servicio que usa drones multirotor para realizar la entrega comercial de paquetes [6]. Así mismo, DHL también inició su proyecto Parcelcopter en 2013, este ha transportado medicamentos a la isla de Juist en el mar del norte [7]. En 2014 Google reveló el proyecto “Wing”, en el que se producen drones con mayor capacidad de carga que los anteriormente mencionados [8].

A partir de estas motivaciones, surgen el FSTSP (Flying Sidekick Traveling Salesman Problem) [9] y posteriormente el VRPD (Vehicle Routing Problem with Drones) como ramas de investigación para este problema. Por estos motivos, es fundamental plantear la posibilidad de usar flotas mixtas, en este caso una de camiones y Vehículos Aéreos No Tripulados, también conocidos como UAVs. Esta investigación pretende contribuir en el estudio de los problemas de ruteo de vehículos con drones, enfocada en la distribución de última milla, un tema relativamente reciente y que cuenta con gran interés por muchas empresas como Amazon [10], DHL [7] y Google [8]. Las contribuciones importantes del presente documento son: la consideración de la incidencia del peso del paquete que lleva el dron, la repercusión que tiene la velocidad del viento en la velocidad del dron; una restricción relevante, la cual es que la batería del dron no exceda el 80% de su capacidad, dado que, si lo hace, este puede tender a caídas. Además, este trabajo evalúa y compara dos soluciones iniciales para el algoritmo ALNS. Estas consideraciones son tomadas en cuenta para lograr presentar una formulación matemática que se aproxime más a la realidad y al mismo tiempo definir el desempeño de los métodos de soluciones iniciales.

Este documento está organizado de la siguiente forma: revisión de la literatura relacionada con VRP y sus variantes (incluido el VRPD) secc. 2.1. En la secc. 2.2, se indica la formulación del modelo matemático, descripción del problema, y la vinculación entre el modelo VRPD y distribución de última milla. Secc. 2.3, diseño del algoritmo de solución, incluyendo dos soluciones iniciales (CWA y CWA Ext.) y algoritmos de destrucción e inserción. Luego, se presenta en la Secc. 3, la experimentación y resultados computacionales, con experimentación de instancia y GAPS del ALNS y resultados característicos para la resistencia de vuelo del UAV y las mejores soluciones para el problema. Finalmente, las conclusiones y recomendaciones para trabajos futuros se proporcionan en las secciones 4 y 5.

2. Metodología

2.1. Revisión de la literatura

La primera investigación referente al VRP se dio con “The Truck Dispatching Problem”, donde se realizó una formulación del problema para la distribución de combustible. En dicho artículo, Dantzig y Ramser [11] obtuvieron una ruta óptima entre una terminal granelera y las terminales de gasolina que atiende esta, mediante el uso de una flota homogénea de camiones repartidores de gasolina. Se estableció la demanda de los productos a entregar y se asignaron las estaciones de abastecimiento de los camiones repartidores, de manera que fuera satisfecha la demanda y que, al mismo tiempo, el recorrido transitado por la flota de vehículos fuera el mínimo posible.

Con el paso de los años se han estudiado distintas variantes a este problema, con el fin de aproximarlos más a la realidad, una de ellas es el uso de una flota heterogénea denominado HVRP [12], donde se analiza el comportamiento de la respuesta al tener vehículos de diferentes capacidades. En este el problema, The Fleet Size and Mix Vehicle Routing Problem[13], realizó una comparación entre distintos algoritmos y finalmente se estableció la combinación óptima en la que se minimizaran los costos totales. En los últimos años, y motivado en gran manera por los cambios climáticos, se han incorporado a la flota de camiones, vehículos con combustión alterna como biodiésel, gas natural e híbridos [14], esto con el propósito de investigar las posibilidades de contribuir a mejorar el medio ambiente y reducir la huella de carbono de las rutas actuales. Luego de esta contribución surge el Problema de Ruteo de Vehículos Eléctricos (EVRP) [15], donde se consideran restricciones como el tiempo de carga de la batería del vehículo y la disponibilidad de la infraestructura de carga a lo largo de la ruta. Otra de las posibilidades para reducir la huella de carbono y además minimizar los tiempos de ruta, es el uso de camiones en combinación con drones para la entrega de paquetes en la distribución de última milla, este problema es también una variante del VRP; se denomina Problema de Ruteo de Vehículos con Drones (VRPD). Murray y Chu en 2015 realizaron un primer acercamiento al problema abordando el denominado FSTSP, el problema consistía en determinar las asignaciones óptimas de clientes para un UAV que trabaja en conjunto con un único camión de reparto; el objetivo del modelo es reducir el tiempo de duración de la ruta de los camiones y drones para dar servicio a todos los clientes y devolver ambos vehículos al depósito. En este problema se contaba con un dron que podía desplazarse sin necesidad de que el camión esperara, el problema fue resuelto mediante un modelo de Programación Lineal Entera Mixta (MILP) y probado con instancias de la literatura [9].

En 2017 [16], presentaron un problema multiobjetivo de ruteo de vehículos con flota heterogénea de UAV, donde se consideró la autonomía limitada de los vehículos, múltiples estaciones de carga y se respetaron los requisitos operativos, se formuló usando un modelo MILP, resuelto mediante un algoritmo de búsqueda dinámica y validado mediante un caso de estudio de dos capas.

A mediados de 2018, los autores [17] estudiaron sobre el uso de los drones para la distribución de última milla por considerarlos una idea prometedora para reducir el tráfico excesivo por carretera. Con el fin de superar las dificultades causadas por la corta capacidad relativa de vuelo de los drones, estos vehículos aéreos no tripulados propusieron un concepto innovador al usar los camiones como plataformas móviles de despegue y aterrizaje, partiendo de que la ruta de los camiones ya estuviese dada, el objetivo es optimizar el cronograma de los drones lanzados desde el camión hasta el préstamo del servicio a un conjunto determinado de clientes, resuelto con dos modelos de programación entero y mixto, probado con la generación de tres tipos de instancias. Durante 2019 los autores [18] consideraron la utilización de los drones con camiones para la distribución de paquetes, teniendo dentro de su modelo restricciones de capacidad, carga y tiempo de vuelo de los drones, además los clientes de este modelo se agrupaban. El objetivo de esta investigación fue minimizar el costo de las rutas, para lograrlo emplearon el algoritmo ALNS (Adaptive Large Neighborhood Search) del cual se obtuvo buenas soluciones, que posteriormente fueron comparadas y evaluadas mediante 112 instancias de [19]. Es por esto por lo que, a partir de las combinaciones y de las restricciones que se han establecido en el área de investigación, se ha logrado abarcar un problema de ruteo de vehículos con drones donde se consideró que la flota de drones es homogénea, la resistencia de la batería dron, su velocidad, y el tiempo de vuelo de este el cual se relaciona directamente con los factores anteriormente mencionados.

2.2. Formulación del modelo matemático

Considerando los hallazgos de la revisión de la literatura específica del área, se plantea un modelo matemático adaptado de aquel formulado en [18], el modelo del presente documento tiene en consideracion aspectos importantes del dron como lo son las implicaciones del peso del paquete, la incidencia del viento, la limitacion del consumo de bateria.

2.2.1.Descripción del problema

El problema de distribución propuesto consta de dos subproblemas a resolver sugerido por [9]. El primero consistió en la necesidad de trazar una ruta que optimizara el uso de una flota heterogénea que llevara a cabo la repartición de bienes y servicios encaminados a la distribución de última milla; el segundo, en determinar cuál ruta reduce los costos operacionales y cumple con restricciones como el tiempo máximo de la ruta.

Para el primer caso, en el trazado de la ruta, es necesario considerar las particularidades de una flota que está compuesta por un conjunto de vehículos y UAVs, entre ellas se tiene que el camión de reparto funciona bajo el mando de un conductor y debe estar equipado por un dron, ambos con la misma misión de entregar paquetes a los clientes, los cuales son servidos una sola vez, ya sea por un camión o por un dron, además, cada camión transporta su correspondiente dron y este debe partir y regresar a un único depósito, el retorno puede ser en tándem o de manera independiente. En los momentos en los que el dron no realiza ninguna entrega, permanecerá en el camión para ahorrar energía de la batería, en cambio, cuando el dron está en operación es lanzado por el vehículo desde un nodo 𝑖 y recogido por este en varios nodos 𝑘, este proceso se realiza por cada uno de los arcos de la ruta trazada. En el modelo se considera que el dron únicamente puede atender un cliente a la vez, debido a que la capacidad de carga útil del dron no es ilimitada y que el peso de determinados paquetes solamente puede ser transportado por los vehículos, puesto que los drones no cuentan con la fuerza idónea para llevar productos que excedan su capacidad de elevación. Así mismo, se tiene un tiempo asociado al lanzamiento y recuperación del dron; igualmente, se considera un tiempo de servicio para entregar los paquetes a los clientes. Por otra parte, los vehículos tienen una restricción vinculada, a fin de que no se exceda un límite de tiempo de operación por día.

2.2.2. Modelo matemático

Índices

𝑖 = Nodo de partida, 𝑖 = 1, 2, 3 … , 𝑛

𝑗 = Nodo de llegada, 𝑗 = 1, 2, 3 … , 𝑛

𝑘 = Nodo de recuperación, 𝑘 = 1, 2, 3 … , 𝑛

𝑛 = Número de nodos

𝑣 = Número de vehículos

𝑑 = Número de drones

𝑠 = Número de tuplas ⟨i, j, k⟩

Listado de conjuntos

𝐶 = Conjunto de clientes {1, 2, 3 … , 𝑐}

𝐶’ ⊆ 𝐶= Subconjunto de clientes visitados por drones

𝑉 = Conjunto de vehículos {1, 2, 3 … , 𝑣}

𝐷 = Conjunto de drones UAV {1, 2, 3 … , 𝑑}

𝑁 = Conjunto de todos los nodos incluyendo el depósito cuando sea 0 {0, 1, 2, … 𝑐, (𝑐 + 1)}

𝑁₀ = Conjunto de nodos desde el cual un vehículo puede partir {1, 2, 3 … , 𝑐}

𝑁_𝑖^∗ = Conjunto de nodos al cual un vehículo puede llegar {1, 2, 3 … , 𝑐, (𝑐 + 1)}

∆⁺(𝑖) = El conjunto de nodos que pueden ser alcanzados desde el nodo 𝑖 ∈ 𝑁0 = 𝑁^∗ {𝑖}

∆⁻(𝑖) = El conjunto de nodos que pueden ser usados para alcanzar el nodo 𝑖 ∈ 𝑁_𝑖^∗ = 𝑁₀ {𝑖}

𝐴 = Conjunto de todos los arcos {(𝑖, 𝑗): 𝑖 ∈ 𝑁0 , 𝑗 ∈ ∆⁺ (𝑖)}

Parámetros determinísticos

𝑞 = Capacidad del vehículo [𝐾𝑔]

𝑑_𝑖 = Demanda del cliente 𝑖 ∈ 𝐶 [𝐾𝑔]

𝑊 = Capacidad de carga del dron [𝐾𝑔]

𝑒 = Capacidad de la batería dada en tiempo de vuelo del dron [𝐻]

𝐴_𝑖𝑗 = Distancia requerida para ir desde el nodo 𝑖 ∈

𝑁0 hasta el nodo 𝑗 ∈ 𝑁_𝑖^∗[𝐾𝑚]

𝑡_i𝑗^𝑣 = Tiempo requerido por un vehículo para viajar desde 𝑖 ∈ 𝑁₀ 𝑎 𝑗 ∈ 𝑁_i^∗ [𝐻]

𝑡_i𝑗^d = Tiempo requerido por un UAV para viajar desde 𝑖 ∈ 𝑁₀ 𝑎 𝑗 ∈ 𝑁_i^∗[𝐻]

𝐶_𝑖𝑗^𝑣 = Costo de operar un vehículo para viajar desde 𝑖 ∈ 𝑁₀ 𝑎 𝑗 ∈ 𝑁_i^∗[$]

𝐶_𝑖𝑗^𝑑 = Costo de operar un UAV para viajar desde 𝑖 ∈ 𝑁₀ 𝑎 𝑗 ∈ 𝑁_i^∗[$]

𝑀 = Un número suficientemente grande

𝑆_𝑒^𝑖𝑣= Tiempo de servicio del vehículo al cliente 𝑖∈𝐶 [𝐻]

𝑆_𝑒^𝑖𝑑= Tiempo de servicio para el UAV al cliente 𝑖∈𝐶’ [𝐻]

𝑆𝐿= Tiempo requerido para lanzar el UAV [𝐻]

𝑆𝑅= Tiempo requerido para recuperar el UAV [𝐻]

𝑇_𝑚𝑎𝑥= Tiempo máximo de duración de una ruta [𝐻]

𝑇_{𝑚𝑎𝑥𝑑}=Tiempo máximo de duración de una ruta para el dron [𝐻]

El 𝐶_𝑖^𝑗𝑣 se definió como el producto entre la distancia recorrida, el consumo y el precio de combustible. El costo del dron se estableció como un factor ∝ del costo del vehículo. Adicionalmente, fue necesaria una nomenclatura complementaria para identificar los posibles puntos de partida o salida de un dron; una salida del dron se compone por tres nodos de la ruta desde los cuales se puede operar. Siendo 𝑃 el conjunto de posibles salidas, representadas por las tuplas ⟨𝑖,𝑗,𝑘⟩. En donde el primer nodo (𝑖) representa la posición de lanzamiento del UAV, el segundo nodo (𝑗) representa al cliente que es visitado por el UAV y, finalmente, el tercer nodo (𝑘) representa la posición de recuperación del UAV por parte del vehículo. Un elemento ⟨𝑖,𝑗,𝑘⟩, pertenece al conjunto 𝑃 si se obedecen las siguientes condiciones:

• •Posición del nodo de partida del dron o lanzamiento 𝑖∈𝑁₀ de una tupla: es la ubicación desde la cual se puede lanzar un UAV y puede salir un camión.

• Posición del nodo llegada del dron/entrega 𝑗∈{𝐶′:𝑗≠𝑖} de una tupla: es el conjunto de clientes que pueden ser atendidos por el UAV diferente de la posición de lanzamiento 𝑖.

• Posición del nodo de encuentro camión-dron o de recuperación 𝑘∈{𝑁𝑖∗∶𝑘≠𝑖,𝑘≠𝑗,𝑆𝐿+𝑆𝑅+𝜏 𝑖𝑗𝐷 + 𝜏 𝑗𝑘𝐷+𝑆𝑒 𝑗𝐷 ≤𝑒} de una tupla: es la ubicación en la que el camión puede recuperar el UAV respetando la duración de la batería.

Conjuntamente, las tuplas ⟨0,𝑖,(𝑐 +1)⟩ se encuentran excluidas del conjunto 𝑃 para todo 𝑖∈𝐶’. Dichas tuplas hacen referencia a las entregas que deben realizar los drones, donde sus posiciones de lanzamiento y recuperación son en el depósito. Para una salida 𝑠=⟨ 𝑖,𝑗,𝑘⟩, se define el costo como 𝐶𝑠𝐷=𝐶𝑖𝑗𝐷+𝐶𝑗𝑘𝐷. Adicionalmente, se establece a 𝑃𝑖+ como todas las salidas de 𝑃 lanzadas desde el nodo 𝑖∈𝑁0, 𝑃_𝑘⁻ como todas las salidas de 𝑃 con un nodo de recuperación 𝑘∈ 𝑁_𝑖^∗ y 𝑃_𝑗 como todas las salidas de 𝑃 que sirven al cliente 𝑗∈𝐶’.

Variables

𝑝 _𝑖^𝑗𝑣: Indica si un cliente 𝑗∈𝐶 se visita después de la ubicación 𝑖∈𝑁𝑜 en la ruta del camión 𝑣∈𝑉 [0 o 1]

𝑥 _𝑖^𝑗𝑣: 1, si el vehículo v∈V viaja desde 𝑖 hasta 𝑗 .

0, en otros casos

𝑦 _𝑆𝑣: 1, si la salida 𝑠∈𝑃 ⟨𝑖,𝑗,𝑘⟩ es usada en la ruta del vehículo 𝑣∈𝑉. 0, en otros casos

𝑧 _𝑖^𝑗𝑣: 1, si la salida 𝑗∈𝐶 es visitada después de la ubicación 𝑖∈𝑁𝑜en la ruta del vehículo 𝑣∈𝑉 0, en otros casos

𝑢 _𝑖^𝑣: Indica la posición de la visita 𝑖∈𝑁 en la ruta del vehículo 𝑣∈𝑉 [(𝑥,𝑦)𝑘𝑚]

𝜏 _𝑖^𝑣: Indica el tiempo de la ruta del vehículo 𝑣∈𝑉 desde un nodo 𝑛 llegando a la ubicación 𝑖∈𝑁 [𝑚𝑖𝑛]

𝜏 _𝑖^𝐷: Indica el tiempo de un dron 𝑑∈𝐷 desde el vehículo 𝑣∈𝑉 llegando a la ubicación 𝑖∈𝑁 [𝑚𝑖𝑛]

Donde:

en [𝑚𝑖𝑛], siendo 𝑉_𝑣 la velocidad del vehículo y 𝑑𝑣 la distancia que recorre el vehículo.

en [𝑚𝑖𝑛], siendo 𝑉𝑑𝑐 la velocidad del dron cuando lleva un paquete o carga, 𝑉_𝑑𝑑 la velocidad del dron cuando no lleva un paquete o carga y 𝑑𝑑 la distancia que recorre el dron.

Función objetivo

a. Costo operacional de los vehículos

b. Costo operacional de los drones UAV

c. F.O. = Min {Costo operacional de los vehículos + Costo operacional de los drones}

La función objetivo (3) minimiza el costo operativo de la visita a los clientes, esta se compone de dos partes: el costo de operación de la flota de vehículos (1) y el costo de operación de la flota de drones (2), utilizados para la distribución de los recursos. Para el caso de los drones se consideran dos costos: el primero, es el costo de operación desde la ubicación 𝑖, en que se puede lanzar un UAV hasta la ubicación de entrega de un cliente 𝑗. El segundo, es el costo de operación desde la ubicación del cliente j hasta el nodo de recolección 𝑘 del UAV. La suma de los dos costos anteriores representa el costo de operación total de los drones y, a su vez, con la suma de los costos (1) y (2) se obtiene el costo total del recorrido de vehículos y drones trazado en la ruta.

Es importante destacar que, a diferencia de lo presentado en el modelo propuesto por [18] donde la velocidad del dron es un valor fijo, en la presente investigación se propone una ecuación para la velocidad, la cual considera el factor de carga del dron y la velocidad promedio del viento. Así mismo, se consultó hasta que porcentaje es recomendable usar la batería del dron, esto con el fin de no forzar el margen de uso de esta. La cual no debe estar por debajo del 30%, ya que puede ocasionar caídas repentinas, deterioro en el funcionamiento del dron o de su batería [20]. Por esta razón se incorpora una restricción sobre la carga energética del dron, que limita su funcionamiento si la capacidad de su batería (dada en tiempo de vuelo), se encuentra por debajo del 30% de la capacidad total. A continuación, se presentan las ecuaciones de velocidad del dron cuando lleva un paquete y cuando no lo está haciendo.

Ecuaciones de velocidad del dron:

Donde:

𝑉_𝑑𝑐 = Velocidad del dron cuando lleva un paquete o carga [𝑚𝑠⁄]

𝑉_𝑑𝑑= Velocidad del dron cuando no lleva un paquete o carga [𝑚𝑠⁄]

𝑚𝑡 = Masa del dron sin batería [𝐾𝑔]

𝑚𝑏= Masa de la bateria del dron [𝐾𝑔]

𝑚𝑙 = Masa de la carga o paquete en el dron [𝐾𝑔]

𝑉_𝑓𝑑 =Velocidad de la ficha del dron [𝑚𝑠⁄]

=Velocidad promedio del viento

[𝑚𝑠⁄]

Las restricciones a las que se encuentra sujetas la función objetivo del problema son las siguientes:

La restricción (6) asegura que cada cliente sea visitado exactamente una vez, ya sea por un vehículo o por un avión no tripulado. La restricción (7) asegura que todos los camiones deben salir del depósito como máximo una vez. De manera similar, en (8) se garantiza que todos los camiones deben regresar al depósito como máximo una vez. Además, está prohibido viajar entre depósitos según lo dado por (9). Las restricciones de conservación de flujo para el vehículo se definen en (10). Las restricciones (11) y (12), se refieren a la eliminación de recorridos secundarios para el vehículo en relación con la posición de un cliente visitado en la ruta del vehículo. También, la restricción de capacidad de los vehículos y drones se da en (13). Igualmente, el UAV puede ser lanzado y recuperado como máximo una vez desde cada nodo, como se indica en las restricciones (14) y (15). La restricción (16) asegura que si el UAV se lanza y recupera en la ubicación 𝑖∈𝑁₀ y 𝑘∈𝑁_∗, el vehículo visita las mismas ubicaciones, respectivamente. La inicialización de los tiempos para el vehículo y el UAV al comienzo de cada ruta se dan en (17) y (18). Así mismo, la duración máxima de una ruta se establece imponiendo un límite al tiempo de regreso al depósito por los vehículos, como se indica en (19) y (20).

Las limitaciones de tiempo para el movimiento del vehículo (21), definen el momento en que el vehículo llega a la ubicación con respecto a las acciones correspondientes que pueden ocurrir en el medio. Las limitaciones de tiempo para el movimiento del dron establecen el tiempo en que un UAV visita un cliente según la posición del vehículo, estas se definen en (22) y (23). La restricción (22) asegura que, si se lanza un UAV desde la ubicación 𝑖∈𝑁0 al cliente 𝑗∈𝐶’, entonces el tiempo de llegada del UAV al cliente 𝑗 debe ser mayor que el tiempo de llegada del vehículo a la ubicación 𝑖 más el tiempo de viaje del UAV entre la ubicación y el tiempo de lanzamiento del UAV. De manera similar para las operaciones de recuperación, la restricción (23) asegura que si un UAV termina de atender a un cliente 𝑗∈𝐶’ y vuela de regreso al vehículo en la ubicación 𝑘∈𝑁∗, la suma de los tiempos de llegada del UAV a la ubicación 𝑗, el tiempo de viaje entre las ubicaciones 𝑗−𝑘 , el de servicio del UAV en el nodo cliente 𝑗 y el tiempo necesario en la recuperación de este; debe ser menor que el tiempo de llegada del UAV a la ubicación 𝑘.

Las restricciones de sincronización de tiempo para el vehículo y el UAV se definen en (24)-(27). Estas restricciones imponen que el lanzamiento y las operaciones de recuperación estén sincronizadas en el tiempo. Teniendo en cuenta que la sincronización no es necesaria cuando el vehículo y el UAV están regresando por separado al depósito, por ende, las restricciones (26) y (28) se definen para el conjunto de clientes 𝐶 en lugar de 𝑁∗. La restricción de resistencia para la batería del UAV viene dada por las variables de lanzamiento y recuperación. El tiempo máximo de vuelo de un dron en una ruta, no puede estar por debajo del 30% de la capacidad de su batería, esto se asegura con la restricción (28). La restricción (29) asegura que si se lanza un UAV desde la posición 𝑖∈𝑁0 al cliente 𝑗∈𝐶 para recuperarlo en la posición 𝑘∈𝑁∗, entonces la diferencia de tiempo entre las operaciones debe satisfacer el tiempo de la duración de la batería del dron.

La variable binaria 𝑝𝑖𝑗𝑣 define el orden en que el vehículo realiza las visitas, estableciendo que un cliente es visitado antes que otro en la ruta. El valor de esta variable se asigna de acuerdo con la diferencia de posiciones entre la ubicación del cliente y el punto de partida del dron, como se indica en (30) y (31). Suponiendo la ubicación de salida 𝑖∈𝑁0 y el cliente 𝑗∈𝐶, si el vehículo visita al cliente 𝑗 después de estar en la ubicación 𝑖, entonces la diferencia entre las variables de posición 𝑢𝑖𝑣 será positiva, imponiendo el valor uno a 𝑝𝑖𝑗𝑣. Por el contrario, si la diferencia es negativa, 𝑝𝑖𝑗𝑣 se impone el valor de cero. Cuando un vehículo no visita a un cliente en la ruta, establecido en las restricciones (11) y (12) se define la posición de la visita como cero. Por lo tanto, en la restricción (31) se agrega un 1 en el lado derecho para garantizar que 𝑝𝑖𝑗𝑣 permanezca en cero cuando tanto 𝑢𝑖𝑣 como 𝑢𝑗𝑣 tomen el valor de cero.

La restricción (32) fue definida para evitar que ocurrieran nuevos lanzamientos mientras el UAV estaba volando en la ruta. La hora de llegada de un UAV a otra ubicación del cliente 𝑏∈𝐶 desde el cual el UAV puede ser lanzado nuevamente tiene que ser mayor que la hora de llegada a la ubicación 𝑘∈ 𝑁∗, sí y solo si un dron se recupera en la ubicación 𝑘. Esto se tiene en cuenta en el orden en que el vehículo visita los lugares 𝑖 y 𝑏 (𝑝𝑖𝑗𝑣), si se produce un lanzamiento desde la posición 𝑏 y si un avión aéreo no tripulado se recupera en la ubicación 𝑘. Finalmente, el dominio de las variables se define en las restricciones (33)-(36).

2.3. Diseño del ALNS

El modelo matemático formulado se resolverá mediante el uso del algoritmo Adaptive Large Neighborhood Search, cuyo procedimiento de solución se presenta en la figura 1. Este algoritmo consiste en mejorar progresivamente una solución inicial, destruyendo y reparando reiterativamente la solución parcial [21].

2.3.1. Algoritmo de generación de la solución inicial

En esta investigación, la solución inicial se dividió en dos pasos, el primero enfocado a la creación de las rutas para los vehículos haciendo uso de la heurística de Clarke y Wright en su forma básica, y de la extensión propuesta por [23] en la que se incluyen los parámetros de 𝜆 y 𝜇; adicionalmente, a esta última se le aplica un operador por desplazamiento (Shift mutación) con el fin de intentar mejorar la solución inicial. En el segundo paso, para la construcción de la solución inicial, se realiza la inserción de los drones a las rutas creadas mediante el algoritmo de inserción de drones.

El algoritmo de ahorros de Clarke y Wright (CWA) [24] se preocupa por el enrutamiento inicial de suministros, los cuales se llevan a cabo desde un depósito central a una gran cantidad de puntos de entrega que, a través de una búsqueda en un espacio delimitado, logra encontrar los ahorros más representativos entre dos puntos. El proceso que sigue el algoritmo se basa en hallar los ahorros de distancia entre cada par de clientes, seguidamente se procede a ordenarlos de mayor a menor, posteriormente se va trazando la ruta entre dos puntos en el sistema, respetando siempre que la suma de las cargas asignadas sea por lo menos igual a la capacidad del vehículo transportador de mercancía.

Se evalúan dos rutas diferentes (0,...,𝑖,0) 𝑦 (0,𝑗,...,0) que pueden ser combinadas formando una nueva ruta (0,...,𝑖,𝑗,...,0), la nueva ruta se crea sí y solo si el ahorro (en distancia) obtenido por dicha unión es mayor que cero. El ahorro para cada conjunto de rutas r es calculado mediante la siguiente expresión:

El algoritmo de ahorros tiene una versión paralela en la que se trabaja sobre todas las rutas simultáneamente, y otra secuencial que construye las rutas una a la vez. En este trabajo de investigación solo se hizo énfasis en la segunda versión.

Como se mencionó previamente, la extensión del algoritmo de ahorros presenta los parámetros de peso 𝜆 y 𝜇. El parámetro 𝜆 es usado para evitar que se generen rutas circulares, impidiendo así dar lugar a la unión de puntos lejanos. [25] y [26], formularon este parámetro entre 0,1≤𝜆≤2. Adicionalmente, se añadió el parámetro 𝜇 el cual tiene en cuenta la asimetría en términos de distancia desde el depósito hasta los nodos que se están evaluando, por ello [23] añadió a la expresión (37) un nuevo termino que estaba entre 0 ≤𝜇≤2. Debido a esto, la ecuación de ahorros sufre unas leves modificaciones, las cuales se presentan a continuación.

2.3.2.Algoritmos de destrucción

Luego de la creación de la ruta inicial donde solo se consideran los vehículos, se procede a aplicar el algoritmo de inserción para incluir a los drones en las rutas creadas. Este algoritmo tiene en cuenta las restricciones de capacidad de vuelo y de carga del dron; el proceso inicia identificando los clientes que pueden ser atendidos por los drones, luego de esto se evalúa en la ruta del vehículo las posibles tuplas que se pueden generar para atender a dichos clientes trazando la que permita reducir los costos totales de la ruta.

El algoritmo sigue su procedimiento aplicando un algoritmo de destrucción a la solución inicial s, el cual consiste en remover un número de nodos de la solución. La cantidad de nodos que se quitan depende de la tasa de destrucción y esta es definida previamente al momento de seleccionar la heurística que va a ser aplicada en cada fase. La tasa es un número entero determinado aleatoriamente entre un intervalo, dicho intervalo se define con un valor mínimo (𝑑𝑚𝑖𝑛) y uno máximo (𝑑𝑚𝑎𝑥), en función del número de clientes totales. La tasa de destrucción debe estar por debajo del 0,5 de n y se sugiere que esta sea entre el 10% y 40% cuando se trate de instancias pequeñas (50 a 100 clientes) y entre 30 y 60 clientes cuando el número de nodos exceda dicho intervalo [27]. Las normas por las cuales se rige la elección de los clientes dependen de las heurísticas empleadas; las heurísticas usadas en la remoción de los clientes fueron, Random Removal (figura 2), Worst Removal (figura 3) y Cluster Removal (figura 4); a continuación, se hace una explicación del funcionamiento de cada una de ellas.

La eliminación aleatoria o Random Removal se encarga de retirar de manera uniforme y aleatoria γ clientes de la solución inicial, dando lugar a rutas parciales.

La heurística de destrucción Worst Removal, Critical Removal o peor eliminación, basa su mecanismo en ordenar los clientes servidos de una solución 𝑠 de manera decreciente de acuerdo con sus costos, estos son guardados en una lista ℒ; aplicando como criterio la diferencia de costos Δ𝐶−𝑖(𝑠) de cada cliente en la ruta, esta se calcula a partir de la diferencia entre las funciones [𝐶(𝑠) -𝐶−𝑖(𝑠)], la cual se compone de los valores de la solución 𝑠 incluyendo y excluyendo al cliente 𝑖. El propósito central de esta heurística es eliminar al cliente de mayor costo. La eliminación en cada iteración se hace mediante el operador ⌊𝜗𝑚⌋ [28], con este se obtiene la posición del cliente a ser retirado de las lista. La forma en la que se realiza esto es eligiendo uno a uno los γ clientes, generando así un número aleatorio 𝜗∈[0,1] y un factor 𝑚≥1 con el que se logra la diversificación de la heurística y controla su aleatoriedad, cuando 𝑚=1 la remoción de los clientes en la solución 𝑠 es hecha totalmente aleatoria.

Finalmente, la heurística de eliminación Cluster Removal inspira sus mecanismos en los resultados observados a partir del comportamiento de la remoción orientada a la distancia [29], para hacerlo divide el conjunto de nodos de una ruta en 𝑘 subconjuntos disjuntos, siendo 𝑘≥2, agrupando los nodos que se encuentren de alguna manera “cercanos entre sí”. El algoritmo Cluster Removal comienza su procedimiento con la tasa de destrucción que es un valor aleatorio entre 10% y 30% del total de clientes, luego se procede a determinar el número de ?? subconjuntos a conformar, este se calcula aproximando hacia abajo la división de los clientes totales y el número de clientes a remover establecidos por la tasa. Los 𝑘 subconjuntos se generan mediante el método de K-means, que tiene en cuenta los centroides con el fin de calcular los clientes que se encuentre más cercanos entre sí y así formar los clústeres, este proceso se repite iterativamente hasta lograr obtener clústeres de forma uniforme en la solución.

Se tienen dos casos: el primero es en el que el número de clientes a remover es mayor que el número de clientes pertenecientes al clúster, para este escenario se remueven todos los clientes que lo componen, la selección de este se hace al azar, luego los clientes faltantes para cumplir la tasa de destrucción son elegidos al azar de otro subconjunto de 𝑘 clústeres creados. Por otro lado, tenemos el caso en que la cantidad de nodos a remover es menor al número de clientes del subconjunto seleccionado, estos son removidos en su totalidad. Estos procesos se repiten hasta que se completa con un número de clientes definidos por la tasa de destrucción. Una ventaja de este método es que contempla distintos clientes de distintas rutas y logra agrupar nodos de instancias de gran tamaño, debido a que el recorrido empieza y finaliza en la posición central y visita varios grupos, estos subconjuntos se implementan para mejorar notablemente la distancia 𝑑 entre las rutas creadas y de esta manera mejorar los costos de solución.

La heurística de destrucción trata de eliminar un grupo completo y asignarlo a un recorrido distinto, por ello, lo que se busca es evitar una larga distancia entre la unión de ellos [30].

2.3.3. Algoritmos de inserción

Para la fase de reparación de las rutas, se proponen 3 algoritmos, de los cuales dos son de tipo voraz, Greedy Insertion with Noise function (GIN) y Regret-n Insertion (nRI) debido a que reparan la solución parcial tomando a los clientes removidos e insertándolos mediante el cálculo del costo de inserción de los clientes que aportan una ganancia a la solución parcial. La tercera heurística usada se denomina Repair Close Insertion la cual contempla las distancias para realizar la inserción toma al nodo más cercano que pertenece a la solución parcial actual.

La primera heurística se conoce como inserción codiciosa con función de ruido o Greedy Insertion with Noise function (GIN) (figura 5), esta es una versión extendida del Greedy Insertion (GI), la cual se diferencia de esta última porque contempla un grado de libertad a la hora de elegir el mejor lugar en el cual se va a posicionar un nodo, con el que se obtiene un nuevo costo, este se utiliza para localizar la ubicación en la cual se va a reinsertar el nodo. Esta extensión incorpora 𝑑̅, 𝛼 𝑦 ε al costo actual de la ruta, en donde 𝑑̅ equivale a la distancia más larga que recorre el camión entre los nodos, si considera insertar a un nodo entrega del camión. De la misma manera, si se contempla la opción de insertar un cliente como nodo de dron, 𝑑̅ será la distancia máxima que recorre un dron. Asimismo 𝛼 representa un parámetro de ruido que es empleado para variar la selección de los nodos a escoger, este toma el valor de 0,1[32]. Adicionalmente, se tiene ε el cual es un número aleatorio entre [-1,1].

El costo actual se compone por los costos de los camiones y drones. Siendo así como el GIN halla la posición y logra reinsertar un nodo con un mínimo costo, el cálculo mediante el cual se obtiene dicho costo es el siguiente:

La segunda heurística usada es Regret-n Insertion (nRI) (figura 6), al igual que la anterior basa su mecanismo en los costos; sin embargo, en esta heurística no solo se tiene en cuenta la mejor posición de inserción, sino que se evalúan las n mejores posiciones para cada cliente removido buscando con esto evitar los problemas presentados en las heurísticas netamente voraces, uno de ellos es que del conjunto de clientes removidos comienzan a ser insertados en la ruta los clientes con mayores ahorros, pero al evaluar los últimos nodos pueden resultar difíciles de insertar llegando a generar soluciones no factibles. La ecuación (40) empleada en esta heurística calcula la diferencia de valor entre la n-ésima y la mejor posición del cliente, gracias a esto en las primeras iteraciones se puede identificar aquellos clientes críticos y los va insertando en las rutas.

La metodología de la tercera heurística de reparación Repair Closest Insertion (RCI), trata de insertar a los nodos removidos mediante el cálculo de distancias, seleccionando el nodo más cercano que pertenece a la solución parcial actual. El procedimiento comienza seleccionando aleatoriamente a un cliente de la lista ℒ de clientes removidos, luego calcula la mejor inserción de este mediante el criterio del nodo más cercano que pertenece a la ruta parcial, siempre y cuando se cumplan las restricciones de capacidad del vehículo y duración de la ruta, además que el nodo de la ruta parcial no tenga más de dos vecinos en la ruta. Dentro de la evaluación también determina si es menos costosa la inserción de un dron, el nodo insertado es removido de la lista ℒ y se continúa con la selección aleatoria de otro nodo. En el caso de que el nodo seleccionado de la lista ℒ no pueda ser insertado, este nodo es agregado al conjunto 𝐷_𝑁 de clientes faltantes por insertar, los cuales se agregan a la solución parcial mediante el algoritmo de Greedy Insertion with Noise Function (GIN).

2.3.4. Parametros del Algoritmos

Para el adecuado funcionamiento del algoritmo ALNS fue necesario definir cinco parámetros importantes, los cuales son: el periodo de actualización de los pesos, los pesos a asignarse a cada heurística tanto de destrucción como de reparación, basados en el éxito de cada heurística con relación a la función objetivo, el intervalo de la tasa de destrucción usado en estas heurísticas, el criterio de aceptación, el cual admite una solución factible peor que la actual y finalmente un factor que disminuye la temperatura en cada iteración y funciona como criterio de parada en el algoritmo.

El primer parámetro es el periodo de actualización de los pesos, según [28] la actualización de los pesos en las heurísticas debe realizarse cada cierto número de iteraciones 𝑝_𝑢, esto se conoce como periodo de actualización o segmento 𝑢.

Los resultados del trabajo investigativo recomiendan ajustar razonablemente dicho periodo con base en la cantidad de heurísticas a usar en el ALNS, ya que si este es un valor pequeño, puede suceder que algunas de las heurísticas usadas no logren ejecutarse desde el principio, esto ocasiona que sus pesos se reduzcan, dicha reducción puede deberse a la falta de pruebas de la heurística y no a un mal desempeño de esta; con el fin de evitar estos problemas y además reducir el esfuerzo computacional que resultaría cuando se produce una modificación, Ropke y Pisinger [31] sugieren que el reajuste de los pesos suceda cada 100 iteraciones cuando se tengan más de diez heurísticas en total [28].

El segundo parámetro son los pesos adaptativos de las heurísticas, que varían según la calidad de la solución alternativa generada, en cada iteración los pesos son modificados [33]. Como se mencionó en el párrafo anterior, los pesos de las heurísticas son actualizados en un periodo 𝑝𝑢 de iteraciones, estos pesos fueron definidos como 𝑤(𝜃_??) 𝑦 𝑤(𝜃_𝑟), donde 𝜃_𝑑{𝑑=1,2,3} es el conjunto de heurísticas de destrucción y 𝜃𝑟{𝑟=1,2,3} el conjunto de heurísticas de reparación. La elección de la heurística se realiza usando The roulette wheel selection principle, donde en cada iteración la probabilidad de seleccionar una heurística se define como:

Inicialmente, la probabilidad de selección es igual para todas las heurísticas, esta se va ajustando en cada periodo de actualización mediante un factor de reacción 0≤𝜌≤1, la puntuación Ψ que depende únicamente del tipo de solución obtenida con las heurísticas usadas, al final de cada iteración esta se va sumando en 𝜋𝑑 para cada heurística de destrucción y en 𝜋_𝑟 en las de reparación, las cuales inicialmente tienen valor cero. Las puntuaciones de Ψ a asignar se definen en la tabla 1, y finalmente 𝜃𝑑 𝑦 𝜃𝑟 cuentan el número de veces que ha sido usada la heurística. Para hacer la actualización de los pesos se usa para las heurísticas de destrucción y reparación respectivamente:

El tercer parámetro del algoritmo es el intervalo de la tasa de destrucción de estas heurísticas, esta se establece que debe estar por debajo del 0,5 de n, se sugiere que esta sea entre el 10% y 40% cuando se trate de instancias pequeñas (50 a 100 clientes) y entre 30 y 60 clientes cuando el número de nodos exceda el intervalo de clientes para instancias pequeñas [34], esto con la finalidad de ampliar el espacio de búsqueda. Sin embargo, si se fija un valor demasiado alto se pueden presentar inconvenientes en la fase de reparación, debido a que las heurísticas de inserción no están diseñadas para reparar completamente una solución.

El cuarto parámetro es el criterio de aceptación, la función de este parámetro es controlar y aceptar las soluciones creadas en cada iteración del algoritmo; para definirlo se usó la heurística de Recocido Simulado [35], [36], en la que se utiliza un parámetro de temperatura 𝑇 que regula la probabilidad de aceptación. Si la solución 𝑆𝑁 generada por las fases de destrucción y reparación da como resultado una solución con menor valor objetivo que la solución actual 𝑆_𝐶, entonces se acepta. Sin embargo, si la solución 𝑆_𝑁 tiene un valor objetivo más alto, se hace uso del algoritmo de recocido simulado, con el que se calcula una probabilidad de aceptación 𝑃𝑟𝑜𝑏(𝑆_𝑁), la cual es comparada con un número generado aleatoriamente entre cero y uno; la solución únicamente es aceptada si el valor de dicha probabilidad es mayor al número aleatorio 𝑅𝑎𝑛𝑑𝑜𝑚(0,1)<𝑃𝑟𝑜𝑏(𝑆𝑁). El cálculo de la probabilidad se realiza con la diferencia entre el costo de las dos soluciones a razón de una temperatura que disminuye dependiendo del tiempo máximo y el tiempo empleado iterativamente, a continuación, se presenta la fórmula para el cálculo de la probabilidad:

Donde 𝑓(𝑆_𝐶) y 𝑓(𝑆_𝑁) son los valores de la función objetivo evaluados con 𝑆_𝐶 y 𝑆_𝑁 respectivamente, además es necesario definir una temperatura inicial 𝑇𝑖𝑛 para 𝑇 que se reduzca linealmente hasta hacerse cero; esto con la finalidad de que también se establezca como un criterio de parada del algoritmo, para esto se controla la temperatura usando el tiempo transcurrido 𝑡𝑝𝑎𝑠𝑠 desde el inicio del algoritmo y se toma 𝑡𝑚𝑎𝑥 como el límite de tiempo impuesto en el algoritmo. Para calcular la nueva temperatura se emplea la siguiente fórmula:

3. Experimentación y resultados

La experimentación y el análisis de los resultados se presentan en esta sección, se estudiará el rendimiento del algoritmo para el VRPD empleado en esta investigación.

3.1. Validación del algoritmo

Varias de las instancias presentadas en distintos artículos no pueden ser usadas en el problema de ruteo de vehículos con drones, debido a que este contempla parámetros como la capacidad de carga del dron, los tiempos asociados al dron, la capacidad de las rutas, el tiempo de vuelo de la batería del dron, entre otros, esto le da un trato especial a este problema. Por lo tanto, evaluar dicho problema con instancias que no contemplen este tipo de parámetros lleva a un escenario que no refleja correctamente las propiedades reales del problema.

Esta investigación evalúa el desempeño del algoritmo ALNS planteado mediante las instancias propuestas por [37]. Estas instancias presentan una distribución uniforme 𝑈(−𝑑,𝑑), las instancias son generadas al azar y cuentan con una acomodación en la que solo se contempla un depósito central ubicado en las coordenadas (0,0). Los clientes se ubican en una cuadrícula de dimensiones 2𝑑 por 2𝑑 alrededor del depósito. Cada instancia es denominada 𝑛.𝑚.𝑡, donde 𝑛 es el número de clientes para el escenario, 𝑚 es dimensión de la cuadrícula y 𝑡 es el nombre genérico del escenario. Cada instancia se encuentra compuesta por tres columnas: coordenada en X, coordenada en Y, y demanda de cada cliente. Por último, se tiene un archivo denominado RouteData, que proporciona la información de los parámetros del modelo usados para la ejecución del algoritmo.

Las instancias planteadas por los autores tienen ocho grandes grupos, los cuales son de distintos tamaños 6, 10, 12, 20, 50, 100, 150 y 200 clientes respectivamente. cada uno de estos sets contiene tres subgrupos de instancias con configuración diferente, teniendo así un total de 112 instancias donde se plantean diferentes escenarios para el problema de ruteo de vehículos con drones. La dimensión de las cuadrículas se establece en gran medida por los drones, dado que estos pueden operar un recorrido de ida y vuelta con un radio de alrededor de 10 millas desde el centro de distribución [6]; por esta razón disponer de un tamaño por encima de ese podría causar implicaciones y escenarios en los que el dron no pueda funcionar de manera independiente, por ende, se obliga a que el camión interactúe con el dron.

Las instancias usadas tienen seis tamaños de cuadrícula, 5, 10, 20, 30 y 40 millas respectivamente; se distribuyen 6,10,12 y 20 clientes en los tamaños de cuadrícula 5,10 y 20, denominadas instancias pequeñas, se tienen cuatro configuraciones para cada número de clientes y cuadrícula, siendo 48 problemas en total. Así mismo, se cuentan con 50, 100, 150 y 200 clientes distribuidos en cuadrículas de 10, 20, 30 y 40 con cuadros de diferentes configuraciones, siendo 64, consideradas instancias grandes; obteniendo así un total de 112 instancias para realizar la validación del algoritmo. Inicialmente, se hicieron pruebas con instancias de 100 y 200 clientes, encontrando que el tiempo computacional era alto; por lo tanto, por temas de eficiencia computacional, se decide tomar una instancia de cada configuración, es decir, 5, 10 y 20 para instancias pequeñas, 10, 20, 30 y 40 para instancias medianas y grandes, probando el algoritmo en 28 instancias en total. Por último, en el archivo RouteData, se proporciona la información de los parámetros del modelo, sus respectivos valores se encuentran en la tabla 2. Cada instancia se corrió 10 veces, el algoritmo fue realizado y ejecutado en un equipo con procesador AMD A9 - 9410 RADEON R5, 5 COMPUTE CORES 2C + 3G y 8 GB de RAM instalada.

3.2 Resultados validación del algoritmo ALNS

Para la validación del algoritmo, además de los parámetros del modelo y de las instancias, fue necesario establecer los parámetros propios del ALNS, para esto se emplearon los presentados por [18], con excepción del valor de 𝜎₄ el cual tiene un valor de 0, pero en esta investigación se propone usar un valor de 1 con el fin de favorecer la diversificación en la búsqueda de nuevas soluciones.

En la tabla 3 se representan los símbolos, parámetros y valores usados para correr las soluciones del recocido simulado, heurísticas de destrucción e iteraciones que hacen parte del algoritmo ALNS, los valores de dichos parámetros son los recomendados por autores como D. Sacramento [18], Ropke y Pisinger [31] que en sus investigaciones encontraron valores que funcionan para encontrar soluciones aceptables, además de reducir el tiempo de ejecución.

3.3 Experimentación de Instancias

El criterio de evaluación para las instancias será el mejor costo global para hacer la distribución a los clientes. En la tabla 4 se presentan los resultados de la experimentación para 28 instancias, en esta se incluyen los mejores valores de 𝑀𝑖𝑛 (𝐶) (función objetivo) obtenidos por el ALNS, las soluciones iniciales obtenidas mediante el Clarke y Wright (CWA), CWA extendido, la mejora porcentual del ALNS se calcula respecto a las soluciones iniciales dadas por el CWA y el CWA extendido. El CWA genera soluciones más aptas que su versión extendida, de manera que el ALNS no consigue obtener mejoras significativas sobre estas soluciones iniciales, esto observa en los porcentajes de mejoras alcanzados. En instancias pequeñas el porcentaje máximo se aproxima al 22% mientras que instancias grandes ronda el 12% para las soluciones del CWA.

En la figura 7 se utiliza el número de veces en que las soluciones aceptadas proporcionaron una nueva mejor solución global, y las veces que se aceptaron las distintas heurísticas de reparación durante la ejecución del algoritmo.

El análisis del criterio de aceptación se hace mediante el recocido simulado, por esto se grafican las soluciones que son aceptadas en una prueba con 1000 iteraciones para las instancias 20.5.1 (figura 7) utilizando los dos métodos de generación de la solución inicial, el CWA y el CWA extendido.

La solución óptima se encuentra en las primeras iteraciones del algoritmo, en dónde, la curva roja representa la mejor solución encontrada dado un número de iteraciones, presentado así un comportamiento plano. Los puntos azules que se registran en el espacio de solución son aquellas soluciones que fueron aceptadas de acuerdo con el criterio del recocido simulado, las cuales tienen costos peores que las mejores actuales y tiene un cambio agresivo, a media que transcurren las iteraciones estos son menos notorios. De acuerdo con la extensión del problema, el mecanismo del recocido simulado acepta propuestas de costos poco favorables al inicio de la búsqueda, con respecto a la mejor solución conocida y a medida que disminuye la temperatura, las variaciones van desapareciendo moderadamente y se llega a un punto en el que solo se aceptan aquellas soluciones ligeramente menores.

3.4 Resultados de los GAP del ALNS

La figura 8, presenta los ahorros de la solución obtenida por la función objetivo y dada a partir del ALNS para el VRPD respecto a la solución inicial generada por Clarke y Wright y su versión extendida para las instancias de 20 clientes. En la mayoría de los casos, se nota cómo la metaheurística ALNS mejora la solución inicial en un 5- 25%, demostrando así la efectividad del algoritmo en relación con que no permanece en los óptimos locales, permitiendo obtener mejores soluciones.

Para instancias pequeñas, el gráfico muestra dispersión respecto al ahorro en la solución inicial. Esto podría darse porque la solución inicial algunas veces halla soluciones casi óptimas, presentando una brecha pequeña en cuanto a la mejor solución generada por el ALNS. Los ahorros bajo el enfoque VRP, es decir, solo camiones, suelen estar entre el 50-70% y, a veces, incluso mayor, esto demuestra que hay una ventaja del uso combinado de camiones y drones, respecto al del VRP.

Los GAP con respecto a la solución inicial (CWA) y bajo el enfoque de solo camiones (VRP) se calculan a partir de la siguiente expresión:

3.5. Resultados característicos

3.5.1. Resultados de vuelo

La figura 9 presenta el porcentaje promedio del uso de la batería del dron, la cual tiende a aumentar a medida que incrementa la dimensión de la cuadrícula que contiene a los clientes, se debe agregar que este aumento está presente en todos los escenarios.

Usualmente, el consumo promedio de la batería del dron no supera el 45%, además se puede observar que para las instancias de 20 clientes no se tienen variaciones significativas entre la solución inicial generada por el Clarke y Wright y el CWA extendido.

3.5.2. Mejores soluciones

En la figura 10 se representan gráficamente las mejores soluciones encontradas por la metaheurística ALNS para el VRPD y VRP en un escenario con 20 y 100 clientes respectivamente. Cabe destacar cómo el escenario final para los camiones se ve afectado por la adición de drones en el problema, logrando así que las rutas de los drones tiendan a acortar significativamente los recorridos del camión, enviando al dron a aquellos clientes difíciles de alcanzar.

3.6. Resultados característicos

En el proceso de entrega otro factor que puede incidir es el tiempo, este se calcula mediante la velocidad del dron y la distancia recorrida por este. El valor que se estableció para la velocidad del dron inicialmente fue de 50 mph y es interesante estudiar qué sucedería si el dron operara a velocidades más bajas a causa de la carga que va a llevar este y la influencia que el aire impone sobre él.

Al disminuir la velocidad del dron, se amplía el tiempo que le va a tomar a este visitar los nodos, considerando que no puede violar las restricciones establecidas para la capacidad de la batería dada en tiempo de vuelo, es posible que no pueda alcanzar a todos los clientes que atiende en las soluciones actuales que abordan el parámetro fijo de velocidad. Por esta razón se decide evaluar la ecuación de velocidad propuesta en esta investigación en la sección 2.2. En la tabla 5 y tabla 6 se tienen las rutas que hace el dron dentro de cada instancia. Siendo 𝑚𝑡 + 𝑚𝑏 la masa del dron (incluyendo batería cargada), 𝑚_𝑙 la masa del paquete, 𝑐_𝑚 el factor de carga, 𝑣_𝑣 ̅la velocidad promedio del viento, 𝑉_𝑓𝑑 la velocidad de la ficha del dron, 𝑑 la distancia empleada en la ruta, 𝑡 el tiempo sin incluir la influencia del viento y el factor de carga y por último el tiempo contemplando estos dos últimos factores.

4. Conclusiones

La influencia que tiene el estudio de los drones en los procesos de distribución de última milla tiene un efecto sumamente positivo. Las nuevas aplicaciones tecnológicas logran desplegar un sin número de intereses en las compañías que a su vez despiertan la atención de los investigadores en plantear tanto cuantiosos modelos matemáticos como estrategias que guíen a la resolución de estos, y, por ende, aportar un amplio beneficio que proviene de la sinergia del VRPD, esto se puede conseguir gracias a la incorporación de los vehículos aéreos no tripulados en las entregas de mercancías, ya que las rutas que se conforman ahorran tiempo, costos y contribuyen a una programación más eficiente.

Es posible evidenciar que la corriente de investigación respecto al uso de drones en la distribución de última milla tiene tres vertientes, la primera consiste en que los camiones lanzan y recuperan a los drones en un mismo nodo; la segunda los camiones lanzan a los drones en un nodo cliente y este es recuperado en otro nodo cliente de la ruta; la tercera es que los drones pueden ser lanzados y recuperados en cualquier punto de la ruta de los vehículos, lo que ofrece una oportunidad para investigar el uso adecuado del dron dentro del ruteo de vehículos.

Debido a que en las primeras fases la metaheurística ALNS cuenta con un algoritmo de ahorros (CWA) el cual logra encontrar buenas soluciones iniciales, esto influye en gran medida a la aplicación de las heurísticas de destrucción y reparación en instancias pequeñas.

De modo que es difícil llegar a mejores soluciones en pocas iteraciones, originando que durante un número significativo de corridas no se encontraran soluciones que fueran aceptadas por tener un costo mayor al generado inicialmente, ocasionando que la temperatura inicial no disminuyera, lo cual generaba que el algoritmo no avanzara pues no se cumplía el criterio de parada de la temperatura final; por esto las instancias pequeñas del problema presentaron serias dificultades para llegar a la temperatura final del recocido simulado y por ende a mejores soluciones.

Se encontró que el CWA extendido encuentra peores soluciones iniciales que su versión original, el algoritmo ALNS ofrece menos soluciones aptas en CWA (independientemente del tamaño de la instancia) esto se debe, en gran media, a que las soluciones iniciales que presenta son como máximo igual de buenas a las del ALNS. Por consiguiente, en pequeñas instancias el máximo porcentaje se acerca al 19,01%, mientras que para instancias grandes se aproxima al 11,4%.

Es posible afirmar que la combinación de las heurísticas de destrucción y reparación en la construcción del algoritmo ALNS favorece a que se explore sin límites el espacio de búsqueda, al contar con un componente de aleatoriedad y la adaptabilidad de estas.

5. Recomendaciones

Abordar en las futuras investigaciones el comportamiento del ALNS en problemas de ruteo de vehículos con drones que impliquen ventanas de tiempo, combinando tanto el servicio orientado al tiempo como el enfoque orientado a la distancia. Utilizando así dentro del algoritmo, heurísticas de remoción, que permitan combinar estos dos enfoques, de tal manera que las heurísticas de reparación logren generar rutas que mantengan un equilibrio entre la cercanía de los clientes y sus periodos de atención.

Investigar sobre posibles algoritmos de búsqueda local simple y operadores de mutación que comparen todas las combinaciones válidas posibles de los mecanismos de intercambio en las rutas, generando así un aporte de mayor rendimiento a las rutas que se le ofrecen al algoritmo ALNS.

Para futuras investigaciones se podría explorar las otras configuraciones que existen entre camiones y drones, planteando además problemas y escenarios más cercanos a la realidad con la consideración de entregas múltiples, entregas parciales y problemas con una flota mixta de drones, identificando si existen cambios al evaluar la función objetivo; además hacer una experimentación de los parámetros del problema, identificando las posibles incidencias de estos en encontrar soluciones factibles y comprobar cómo influyen en el costo.

Realizar un diseño de experimentos para evaluar el efecto que tiene la variación de algunos parámetros del algoritmo como la tasa de destrucción, el factor de reacción y los puntajes asignados cuando una solución es aceptada.

Referencias

Brown Eric, “E-commerce spurs innovation in last- mile logistics”, MIT News, 2019, Consultado: jul. 07, 2022. [En línea]. Disponible en: https://news.mit.edu/2018/mit-e-commerce-spurs-innovations-last-mile-logistics-0904

“Ecommerce: Evolución y perspectivas 2020 | COMUNICAWEB”. [En línea]. Disponible en: https://comunica-web.com/blog/marketing-digital/ecommerce-datos-evolucion/

B. Balcik, B. M. Beamon, K. Smilowitz, “Last mile distribution in humanitarian relief”, J. Intell. Transp. Syst. Technol. Planning, Oper., vol. 12, núm. 2, pp. 51-63, abr. 2008, doi: https://doi.org/10.1080/1547245080202332

B. Minten, B. Koru, D. Stifel, “The last mile(s) in modern input distribution: Pricing, profitability, and adoption”, Agric. Econ., vol. 44, núm. 6, pp. 629-646, nov. 2013, doi: https://doi.org/10.1111/AGEC.12078

“El nuevo mapa de rescates con drones de DJI - ITSitio”. https://www.itsitio.com/co/nuevo-mapa-rescates-drones-dji/

C. Rose, “Jeff Bezos de Amazon mira hacia el futuro - CBS News”, 2013. https://www.cbsnews.com/news/amazons-jeff-bezos-looks-to-the-future

Bonn, “DHL Press Release English”, 2014. https://www.reutersevents.com/supplychain/3pl/dhl-parcelcopter-launches-initial-operations-research-purposes

“Google tests drone deliveries in Project Wing trials - BBC News”. https://www.bbc.com/news/technology-28964260

C. C. Murray, A. G. Chu, “The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery”, Transp. Res. Part C Emerg. Technol., vol. 54, pp. 86-109, 2015, doi: https://doi.org/10.1016/J.TRC.2015.03.005

I. or its affiliates Amazon.com, “Amazon.com: Prime Air”, 2016. https://www.aboutamazon.com/news/transportation/amazon-prime-air-prepares-for-drone-deliveries

G. B. Dantzig, J. H. Ramser, “The Truck Dispatching Problem”, Source Manag. Sci., vol. 6, núm. 1, pp. 80-91, 1959, doi: https://doi.org/10.1287/mnsc.6.1.80

S. Erdoĝan, E. Miller-Hooks, “A Green Vehicle Routing Problem”, Transp. Res. Part E Logist. Transp. Rev., vol. 48, núm. 1, pp. 100-114, 2012, doi: https://doi.org/10.1016/J.TRE.2011.08.001

B. Golden, A. Assad, L. Levy, F. Gheysens, “The fleet size and mix vehicle routing problem”, Comput. Oper. Res., vol. 11, núm. 1, pp. 49-66, 1984, doi: https://doi.org/10.1016/0305-0548(84)90007-8

S. Erdoĝan , E. Miller-Hooks, “A Green Vehicle Routing Problem”, Transp. Res. Part E Logist. Transp. Rev., vol. 48(1), pp. 100-114, 2012, doi: https://doi.org/10.1016/J.TRE.2011.08.001

J. Lin, W. Zhou, O. Wolfson, “Electric Vehicle Routing Problem”, Transp. Res. Procedia, vol. 12, pp. 508-521, ene. 2016, doi: https://doi.org/10.1016/J.TRPRO.2016.02.007

B. N. Coelho et al., “A multi-objective green UAV routing problem”, Comput. Oper. Res., vol. 88, pp. 306-315, dic. 2017, doi: https://doi.org/10.1016/J.COR.2017.04.011

N. Boysen, D. Briskorn, S. Fedtke, S. Schwerdfeger, “Drone delivery from trucks: Drone scheduling for given truck routes”, Networks, 2018, doi: https://doi.org/10.1002/net.21847

D. Sacramento, D. Pisinger, S. Ropke, “An adaptive large neighborhood search metaheuristic for the vehicle routing problem with drones”, Transp. Res. Part C Emerg. Technol., vol. 102, pp. 289-315, 2019, doi: https://doi.org/10.1016/J.TRC.2019.02.018

“Existencias: Vehicle Routing Problem with Drones Instances”. https://zenodo.ups.edu.ec/Record/zen-1568853

“Uso y cuidado de baterías para drones”. https://elvuelodeldrone.com/blog-de-drones/uso-cuidado-baterias-para-drones/

P. Shaw, “Using constraint programming and local search methods to solve vehicle routing problems”, Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics), vol. 1520, pp. 417-431, 1998, doi: https://doi.org/10.1007/3-540-49481-2_30

W. Li, Y. Wu, P. N. R. Kumar, y K. Li, “Engineering Optimization Multi-trip vehicle routing problem with order release time Multi-trip vehicle routing problem with order release time”, Optim. Ing., 2019, doi: https://doi.org/10.1080/0305215X.2019.1642880

H. Paessens, “The savings algorithm for the vehicle routing problem”, Eur. J. Oper. Res., vol. 34, núm. 3, pp. 336-344, 1988, doi: https://doi.org/10.1016/0377-2217(88)90154-3

T. Pichpibul, R. Kawtummachai, “A heuristic approach based on Clarke-Wright algorithm for open vehicle routing problem”, Sci. World J., vol. 2013, 2013, doi: https://doi.org/10.1155/2013/874349

T. J. Gaskell, “Bases for Vehicle Fleet Scheduling”, J. Oper. Res. Soc., vol. 18, núm. 3, pp. 281-295, sep. 1967, doi: https://doi.org/10.1057/jors.1967.44

P. C. Yellow, “A Computational Modification to the Savings Method of Vehicle Scheduling”, J. Oper. Res. Soc., vol. 21, núm. 2, pp. 281-283, 1970, doi: https://doi.org/10.1057/JORS.1970.52

, S. Ropke, “A general heuristic for vehicle routing problems”, Comput. Oper. Res., vol. 34, núm. 8, pp. 2403-2435, 2007, doi: https://doi.org/10.1016/J.COR.2005.09.012

S. Ropke, D. Pisinger, “An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows”, Transp. Sci., vol. 40, núm. 4, pp. 455-472, 2006, doi: https://doi.org/10.1287/trsc.1050.0135

R. Lutz, “Adaptive Large Neighborhood Search”, 2014, [En línea]. Disponible en: http://dx.doi.org/10.18725/OPARU-3237

S. Ropke , “Heuristic and exact algorithms for vehicle routing problems”, 2006. [En línea]. Disponible en: https://backend.orbit.dtu.dk/ws/portalfiles/portal/3155278/Heuristic+and+exact+algorithms+for+vehicle+routing+problems_Ropke.pdf

S. Ropke , D. Pisinger, “An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows”, Transp. Sci., vol. 40, núm. 4, pp. 455-472, 2006, doi: https://doi.org/10.1287/TRSC.1050.0135

A. T. Phan, T. D. Nguyen, Q. D. Pham, “Traveling salesman problem with multiple drones”, ACM Int. Conf. Proceeding Ser., pp. 46-53, 2018, doi: https://doi.org/10.1145/3287921.3287932

M. Keskin, B. Catay, “Partial recharge strategies for the electric vehicle routing problem with time windows”, Transp. Res. Part C Emerg. Technol., vol. 65, pp. 111-127, abr. 2016, doi: https://doi.org/10.1016/j.trc.2016.01.013

D. Pisinger, S. Ropke , “A general heuristic for vehicle routing problems”, Comput. Oper. Res., vol. 34, núm. 8, pp. 2403-2435, 2007, doi: https://doi.org/10.1016/j.cor.2005.09.012

S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, “Optimization by simulated annealing”, Science, vol. 220, núm. 4598, pp. 671-680, may 1983, doi: https://doi.org/10.1126/science.220.4598.671

V. Černý, “Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm”, J. Optim. Theory Appl., vol. 45, núm. 1, pp. 41-51, 1985, doi: https://doi.org/10.1007/BF00940812

D. Sacramento , “Vehicle Routing Problem with Drones Instances”, Transportation Research Part C: Emerging Technologies, ago. 2018, doi: https://doi.org/10.5281/zenodo.2572764

Notas