El problema del m-Anillo Estrella Biobjetivo

Universidad de Sevilla Facultad de Matem´aticas Trabajo Fin de Estudios Grado en Matem´aticas El problema del m-Anillo Estrella Biobjetivo Valentina L´opez Garc´ıa Tutor Miguel ´ Angel Pozo Monta˜no Departamento de Estad´ıstica e Investigaci´on Operativa. Resumen El Problema del m-Anillo Estrella Biobjetivo (biobjective multi-Ring Star Problem - bmRSP) busca el menor recorrido que sirve a todos los nodos de un grafo completo G= (V, A) siendo Vel conjunto de nodos y Ael conjunto de arcos que los conectan. El conjunto Vest´a formado por el dep´osito central, los clientes y los nodos de transici´on (que ayudan a mejorar las conexiones entre nodos). El recorrido estar´a formado por mciclos partiendo del dep´osito central junto con un conjunto de asignaciones que conectan los nodos que no est´an incluidos en los ciclos. Al minimizar el coste de los anillos y el de las asignaciones, llegamos a este problema biobjetivo. Se encuentran aplicaciones del bmRSP en diversos ´ambitos como la log´ıstica, cadenas de distribuci´on, optimizaci´on de rutas y planificaci´on de redes; concretamente, en optimizaci´on de rutas de transporte, redes de telecomunicaciones, sistemas de transporte p´ublico, redes de cadena de suministro y en log´ıstica urbana. Esta memoria aborda las formulaciones de los diferentes problemas que lo preceden y estudia m´etodos de resoluci´on para formulaciones de programaci´on lineal entera mixta del bmRSP. Concluiremos implementando en Python los m´etodos y algoritmos presentados, comparando algunos de ellos y representando las soluciones obtenidas. iii iv Abstract Let G= (V, A) be a complete graph where Vis the set of nodes and Ais the set of arcs that connect them. The set Vincludes the central depot, the customers, and the transition nodes (aid in enhancing connections between nodes). The goal is to find the optimal tour that visits nodes through mcycles starting from the central depot along with a set of assignments connecting the nodes not included in the cycles. By minimizing the ring costs and the assignment costs, we arrive at the biobjective m-Ring Star Problem (bmRSP). This problem finds applications across various domains such as logistics, distribution chains, route optimization and network planning; specifically, in the optimization of transportation routes, telecommunications networks, public transportation systems, supply chain networks, and urban logistics. This work addresses the formulations of the different problems that precede it and studies resolution methods for mixedinteger linear programming formulations of the bmRSP. We will conclude by implementing in Python the methods and algorithms presented, comparing some of them, and representing the solutions obtained. v vi ´ Indice general Resumen iii Abstract v ´ Indice vii Prefacio ix 1 Preliminares 1 1.1 Definici´ondelproblema................................ 1 1.2 Revisi´ondelaliteratura ............................... 2 1.3 Elproblemadelviajante............................... 4 1.4 Formulaci´ondelRSP................................. 5 1.5 Formulaci´ondelCVRP................................ 6 1.6 Formulaci´ondelCmRSP............................... 8 1.6.1 Formulaci´on de dos ´ındices . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.6.2 Formulaci´on como problema de flujo . . . . . . . . . . . . . . . . . . . . 11 1.7 El problema multiobjetivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 M´etodos de resoluci´on 15 2.1 M´etodobiobjetivo .................................. 15 2.1.1 Obtenci´on de los puntos no dominados soportados . . . . . . . . . . . . . 15 2.1.2 Obtenci´on de los puntos no dominado no soportados . . . . . . . . . . . 17 2.2 Algoritmo de ramificaci´on, acotaci´on y corte . . . . . . . . . . . . . . . . . . . . 19 2.2.1 Divideyvencer´as............................... 19 2.2.2 Enumeraci´on impl´ıcita . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.3 Planosdecorte ................................ 24 2.3 Algoritmo heur´ıstico de dos fases . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3 Implemetaci´on 29 3.1 Generaci´ondeinstancias............................... 29 3.2 Carga de datos y representaci´on de soluciones . . . . . . . . . . . . . . . . . . . 30 3.3 Descripci´on general de la resoluci´on . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Algoritmo de ramificaci´on, acotaci´on y corte . . . . . . . . . . . . . . . . . . . . 31 vii 3.5 Algoritmo heur´ıstico de dos fases . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.6 Algoritmo ´epsilon-restricci´on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4 Ap´endice 45 4.1 Formulaciones no dirigidas del RSP y el CmRSP . . . . . . . . . . . . . . . . . . 45 Conclusiones 49 Lista de figuras 51 Lista de tablas 52 Bibliograf´ıa ......................................... 54 viii Prefacio El Problema Biobjetivo del Multi-Anillo Estrella (biobjective m-Ring Star Problem o bmRSP) es la versi´on biobjetivo del Problema del m-Anillo Estrella. ´ Este consiste en dise˜nar una red de ciclos o anillos que pasan por un dep´osito central y satisfacen las necesidad de ciertos clientes. No se impone la necesidad de incluir a todos los clientes en alg´un ciclo. En su lugar, se asignar´an aquellos que no pertenezcan a un ciclo a un cliente que s´ı forme parte de la red de anillos. De esta asignaci´on surge la posibilidad de considerar dos objetivos a minimizar por separado: el coste de ciclo y el de asignaci´on. Resulta de utilidad a˜nadir puntos de transici´on para ahorrar costes de conexi´on. El problema se puede formular con programaci´on lineal entera mixta, tomando un grafo cuyos nodos representan a los clientes y los puntos de transici´on. Se fijan unos costes de anillo y de asignaci´on a cada arco del grafo y se minimiza el coste de los arcos que son seleccionados para formar la red. Esta memoria tiene como objetivos introducir varios problemas que dan lugar al Problema Biobjetivo del Multi-Anillo Estrella, estudiar distintos m´etodos de resoluci´on para problemas de programaci´on lineal entera e implementar en Python los modelos y algoritmos propuestos. Con el fin de lograr los objetivos anteriores, se ha estructurado el trabajo en 3 cap´ıtulos. El primer cap´ıtulo introduce el bmRSP, proporcionando un an´alisis detallado de la evoluci´on de los problemas que lo preceden y formulando algunos de ellos: el Problema del Viajante, el Problema del Anillo Estrella (Ring Star Problem o RSP), el Problema de Ruteo de Veh´ıculos (CVRP), la versi´on uniobjetivo de nuestro problema y el bmRSP. Se presentan dos maneras de formularlo: como problema de flujo o a trav´es de la eliminaci´on de subtours. En el segundo cap´ıtulo se presentan m´etodos de resoluci´on para problemas de programaci´on lineal entera. Se desarrollan tres m´etodos: el algoritmo de la ε-restricci´on para la obtenci´on de los puntos no dominados soportados y no soportados, el algoritmo de ramificaci´on, acotaci´on y corte, y un algoritmo heur´ıstico de dos fases. El algoritmo heur´ıstico se basa en la resoluci´on del problema de ruteo una vez se han establecido las preasignaciones de los nodos. En el tercer cap´ıtulo, se implementa en Python gran parte del contenido de los ix 6 Secci´on 1.5. Formulaci´on del CVRP m´ın X (i,j)∈A cijxij +X (i,j)∈A cijyij (1.3a) s.a. X (i,j)∈A xij =yjj ∀j∈V, (1.3b) X (i,j)∈A xij =yii ∀i∈V, (1.3c) X j∈V yij = 1 ∀i∈V\ {1},(1.3d) yij ≤X j′∈V xjj′∀i, j ∈V\ {1},(1.3e) X j∈V xij ≤ |S| − 1∀S⊂V: 2 ≤ |S|,(1.3f) xij, yij ∈ {0,1} ∀(i, j)∈A. (1.3g) Las dos primeras restricciones, aseguran que para cada arco del ciclo, tanto el nodo de salida como el de entrada est´e asignado a ´el mismo. En (1.3d) obligamos a cada nodo a estar asignado con un s´olo nodo, que puede ser ´el mismo. Cada v´ertice debe estar o bien en el ciclo o bien asignado a un v´ertice del ciclo gracias a la restricci´on (1.3f). Como en el problema anterior, (1.3f) elimina los subtours. Por ´ultimo, la restricci´on (1.3g) proporciona la naturaleza discreta de nuestras variables. Estos dos ´ultimos modelos consideran que el grafo es dirigido, ya que tenemos variables y costes distintos en los arcos (i,j) y (j,i). Existe la formulaci´on de este problema usando grafos no dirigidos. No obstante, este TFG se centra en la versi´on dirigida del bmRSP y, por tanto, la versi´on no dirigida se omite en esta secci´on. Dado su inter´es te´orico, se incluir´a en el Ap´endice. En el caso en el que el coste de asignaci´on fuera dij =cij, la funci´on objetivo ser´ıa: m´ın X (i,j)∈A cijxij +X (i,j)∈A dijyij 1.5 Formulaci´on del CVRP Como extensi´on del TSP, el Problema de Ruteo de Veh´ıculos busca determinar el mejor conjunto de rutas que debe realizar una flota de veh´ıculos para atender a un conjunto determinado de clientes. Una de sus variantes m´as comunes es el Problema de Ruteo de Veh´ıculos Capacitado (CVRP) que introduce veh´ıculos con capacidad de carga limitada y, en ocasiones, restricciones de duraci´on o distancia. Leite (2023) presenta la formulaci´on de flujo de veh´ıculos de tres ´ındices del libro Toth and Vigo (2002) para el CVRP con restricciones de Cap´ıtulo 1. Preliminares 7 carga. Dado un grafo completo G= (V, A) donde V ser´a el conjunto de nodos que incluir´a a los clientes y al dep´osito central por el cual tienen que pasar todos los veh´ıculos y A ser´an los arcos dirigidos del grafo. El conjunto de veh´ıculos Kdar´a lugar a |K|ciclos. Usaremos 3 ´ındices: i y j para referirnos a los nodos y k para los veh´ıculos. En este problema, hay un ´unico nodo dep´osito que se corresponde con el ´ındice 0 y todos los veh´ıculos tienen la misma capacidad Q. Cada cliente itendr´a una demanda/peso qique deber´a satisfacer alg´un veh´ıculo sin que se exceda la capacidad m´axima del veh´ıculo. Consideramos la variable binaria xijk que tomar´a el valor 1 si el veh´ıculo krecorre el arco que conecta el nodo icon el nodo j. La variable binaria yik indica si la demanda del nodo ies cubierta por el veh´ıculo k. Tendremos un costo cij asociado a cada arco y nuestro objetivo ser´a minimizar el coste del recorrido total. As´ı, la formulaci´on del problema ser´a: m´ın X k∈KX i∈VX j∈V cijxijk (1.4a) s.a. X k∈KX j∈V xijk =X k∈KX j∈V xjik = 1 ∀i∈V\ {0},(1.4b) X j∈V xijk =X j∈V xjik =yik ∀i∈V\ {0}k∈K, (1.4c) X i∈V qiyik ≤Q∀k∈K, (1.4d) X k∈KX j∈V x0jk =X k∈KX i∈V xi0k=|K|,(1.4e) X i∈SX j /∈S xijk ≥yhk ∀S⊂V\ {0}, h ∈S, k ∈K, (1.4f) xijk ∈ {0,1} ∀i∈V, j ∈V, k ∈K, (1.4g) yik ∈ {0,1} ∀i∈V, k ∈K. (1.4h) La funci´on objetivo (1.4a) minimiza el coste total de los arcos usados en el recorrido. La restricci´on (1.4b) asegura que cada cliente ies visitado una vez, por lo que existe un ´unico arco que llega a ´el y otro que parte de ´el. La restricci´on (1.4c) impone que si existe alg´un arco recorrido por el veh´ıculo kque entre o salga del nodo i, entonces la demanda qide ese nodo ser´a cubierta por ese veh´ıculo. La restricci´on (1.4d) asegura que la demanda total asignada a cada veh´ıculo no exceda su capacidad. Como queremos hacer uso de todos los veh´ıculos, la restricci´on (1.4e) asegura que se formar´an |K|recorridos que empiezan y terminan en el dep´osito. Por ´ultimo, la restricci´on (1.4f) elimina los subtours que no contengan al dep´osito 8 Secci´on 1.6. Formulaci´on del CmRSP central. La naturaleza binaria de las variables viene dada por las restricciones (1.4g) y (1.4h). 1.6 Formulaci´on del CmRSP Si modificamos el RSP de manera que podamos tener manillos en lugar de uno s´olo, obtenemos el Problema del m-Anillo Estrella. Tambi´en se puede ver como una extensi´on del CVRP, conservando sus k=manillos y la funci´on del dep´osito central pero con la libertad de no tener que incluir a todos los nodos en el recorrido, asignando aquellos que excluimos a otros que pertenezcan a alguno de los anillos. En este problema podemos incluir nodos Steiner: unos puntos de transici´on que nos ayudar´an a ahorrar costes de conexi´on. Tanto los clientes como los nodos Steiner ser´an representados por nodos y los arcos representar´an las conexiones que hay entre ellos. Lana (2016) presenta dos formulaciones distintas del problema que veremos en esta secci´on. Una de ellas se basa en la eliminaci´on de subtours en la cual volveremos a centrarnos en el caso del grafo dirigido, dejando la versi´on no dirigida en el Ap´endice. La segunda formulaci´on se basa en la conservaci´on del flujo y s´olo est´a formulada para grafos mixtos donde los anillos no son dirigidos. Llamaremos nodos visitados a aquellos que formen parte de un ciclo y nodos asignados a v´ertices a los que tengan conexi´on directa con un v´ertice. Los puntos de transici´on ser´an los nodos Steiner que son visitados y pueden tener alg´un nodo asignado (no es necesario incluirlos todos). Una vez m´as, el objetivo de este problema es minimizar la suma de los costes que provocan la construcci´on de anillos y la asignaci´on de nodos. Al igual que en el CVRP, pondremos un l´ımite superior Q al n´umero de clientes que asociaremos a cada anillo, incluyendo visitados y asignados a sus v´ertices. ´ Esta ser´a la capacidad del anillo, que tiene un s´ımil en el l´ımite real de las conexiones de telecomunicaciones, donde un haz de fibras s´olo puede dar conexi´on a un n´umero limitado de clientes. Volveremos a considerar un grafo dirigido G= (V, A). Tomaremos como nodo 0 el dep´osito central y a˜nadiremos un nodo n+ 1 como r´eplica de ´este para facilitarnos la formulaci´on del problema. Los conjuntos UyWrepresentar´an los clientes y los nodos Steiner respectivamente. Por tanto, tendremos V′=U∪Wconjunto de nnodos ya que no est´a incluido el dep´osito central ni su copia. En total, tendremos V={0, n + 1} ∪ V′con |V|=n+ 2. Tendremos el conjunto de posibles conexiones A={(i, j) : i∈U, j ∈Ci}donde Ci⊂V′ representa el conjunto de v´ertices a los que se puede conectar el cliente i. Supondremos que todos los v´ertices tendr´an un arco con el que se asignan a s´ı mismos, es decir, si i∈Ues un v´ertice entonces i∈Ci. Tendremos en cuenta varias definiciones: Cap´ıtulo 1. Preliminares 9  Camino simple entre i1eik: Sucesi´on de v´ertices y aristas tal que de cada uno de los v´ertices parte una arista hacia el v´ertice sucesor, que comienza en i1y termina en ik.  Anillo R: Un par (E’,A’) donde E′⊂Ees un camino simple del nodo 0 al nodo n+ 1 y A′⊂Ason las conexiones directas entre clientes de U\V(E′) (clientes no visitados por R) y v´ertices de V(E′) (clientes visitados por R).  Asociaci´on de clientes a anillos: Un cliente iest´a asociado a un anillo R si es visitado por el camino simple (i∈V(E′)), o bien est´a conectado a un v´ertice del anillo (existe un v´ertice jtal que (i, j)∈A′).  Coste de conexi´on: coste ce=cij no negativo que se le asigna a cada arista e={i, j} ∈ E.  Coste de conexi´on directa: coste dij no negativo que se le asigna a cada arco (i, j)∈A. Para todo i∈U,dii = 0. Como cada anillo tiene una capacidad Q, para que el problema sea factible se debe cumplir que mQ ≥ |U|. En la Figura 1.2 podemos observar la soluci´on factible de una instancia con n= 26, m = 2, Q = 16,|U|= 23 y |W|= 3: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 26 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 0 Funci´on objetivo: 53.1 Dep´osito Central Clientes Nodos Steiner Figura 1.2: Una soluci´on factible para un problema CmRSP con Q=16. 1.6.1 Formulaci´on de dos ´ındices Esta formulaci´on hace uso de las variables xij eyij que hemos visto en anteriores formulaciones. A˜nadiremos para cada j∈W, la variable binaria wjque vale 1 si y solo si 10 Secci´on 1.6. Formulaci´on del CmRSP el nodo Steiner j est´a en un anillo (ya que no es necesario usar todos los nodos de transici´on). Recordemos que V′=U∪WyV={0, n + 1} ∪ V′. La formulaci´on del grafo dirigido ser´a: m´ın X (i,j)∈A cijxij +X (i,j)∈A dijyij (1.5a) s.a. X j∈V x0j=m, (1.5b) X i∈V xi0= 0,(1.5c) X i∈V xi(n+1) =m, (1.5d) X j∈V x(n+1)j= 0,(1.5e) X j∈V xij =yii ∀i∈U, (1.5f) X j∈V xji =yii ∀i∈U, (1.5g) X j∈V xij =wi∀i∈W, (1.5h) X j∈V xji =wi∀i∈W, (1.5i) X j∈V′ yij = 1 ∀i∈U, (1.5j) X i∈SX j∈V\S xij ≥1 QX i∈UX j∈S yij ∀S⊆V′:S=∅(1.5k) X i∈V\SX j∈S xij ≥1 QX i∈UX j∈S yij ∀S⊆V′:S=∅(1.5l) xij ∈ {0,1} ∀(i, j)∈E, (1.5m) yij ∈ {0,1} ∀(i, j)∈A, (1.5n) wi∈ {0,1} ∀i∈W, (1.5o) La expresi´on (1.5a) representa la funci´on objetivo que minimiza el coste total obtenido como la suma de los costes de conexi´on de los arcos de los anillos y los costes de conexi´on directa de los clientes que se conectan a los nodos de los anillos. Las restricciones (1.5b) y (1.5c) aseguran que desde el nodo 0 salgan marcos pero ninguno llegue a ´el. Las restricciones (1.5d) y (1.5e) aseguran que lleguen marcos al nodo n+ 1 pero ninguno parta de ´el. Las restricciones Cap´ıtulo 1. Preliminares 11 (1.5f), (1.5g), (1.5h), (1.5i) imponen para cada nodo i∈V′=U∪Wque pertenece a uno de los anillos, existe un ´unico arco que incide en ´el y otro que sale de ´el. Las restricciones (1.5j) establecen que un cliente io bien est´a en un anillo (yii = 1) o est´a asignado a un nodo de alg´un anillo (yij = 1 para alg´un j∈Ci). Las restricciones (1.5k) y (1.5l) son las restricciones fraccionales de capacidad del anillo. Estas restricciones, junto a la integridad de las variables xey, imponen que, para un subconjunto Sde nodos, al menos Pi∈UPj∈S zij Qanillos son necesarios para visitar a los clientes asignados a los nodos en S, donde ⌈.⌉indica la parte entera del n´umero redondeado hacia arriba. 1.6.2 Formulaci´on como problema de flujo Asociaremos a cada arista (i, j) sus correspondientes variables de flujo fij, fji. Para esta formulaci´on, ambos arcos tienen el mismo coste cij. Una soluci´on se considerar´a factible si el flujo total en cualquier arista del anillo es exactamente Q. Las variables xidentifican cada anillo a trav´es de dos caminos de flujo. En el primer camino, que va del nodo 0 al n+ 1, el flujo fij representa el n´umero de usuarios que son servidos por el anillo tras visitar el nodo i. En el camino opuesto, que va del nodo n+ 1 al nodo 0, la variable fji representa el n´umero potencial de clientes que pueden ser servidos en el camino de 0 a j. As´ı, el flujo total que sale del nodo 0 es |U|, mientras que el flujo que emana del nodo n+ 1 es mQ. Si un nodo en un anillo tiene un n´umero νde unidades de flujo, absorbe νunidades de flujo que vienen del nodo 0 y νunidades de flujo que vienen del nodo n+ 1. 12 Secci´on 1.6. Formulaci´on del CmRSP m´ın 1 QX (i,j)∈E:i<j cij(fij +fji) + X (i,j)∈A dijyij (1.6a) s.a. X j∈V′ f0j=|U|,(1.6b) X j∈V′ f(n+1)j=mQ, (1.6c) X i∈V′ fi(n+1) = 0,(1.6d) X i∈V′ fi0=mQ − |U|,(1.6e) X j∈V (fij +fji) = 2Qyii ∀i∈U, (1.6f) X j∈V (fij +fji) = 2Qwi∀i∈W, (1.6g) X j∈V′ yij = 1 ∀i∈U, (1.6h) X i∈V (fij −fji)=2X i∈U yij ∀j∈V′(1.6i) fij +fji ∈ {0, Q} ∀(i, j)∈E, (1.6j) fij ≥0∀i, j ∈V, (1.6k) yij ∈ {0,1} ∀(i, j)∈A, (1.6l) wi∈ {0,1} ∀i∈W, (1.6m) La expresi´on (1.6a) es la funci´on objetivo que minimiza los costes de los anillos m´as la de conexi´on directa, dependiendo de las variables de flujo. La restricci´on (1.6b) impone que el flujo total que sale del nodo 0 es igual al n´umero total de clientes, mientras que la restricci´on (1.6c) impone que del nodo n+1 el flujo saliente total es igual a la capacidad total de los m anillos. Las restricciones (1.6d) y (1.6e) indican que no entra flujo al nodo n+1 y que el flujo que entra en el nodo 0 es igual a la capacidad residual(no usada) de los anillos, respectivamente. Las restricciones (1.6f) y (1.6g) imponen que la suma del flujo que entra y sale de todos los nodos de los anillos sea 2Q, asegurando que Q unidades de flujo viajan por cada arista de los anillos. Esto se tiene tanto para los clientes del anillo (yii = 1) como para los nodos Steiner usados (wi= 1). Las restricciones (1.6h) aseguran que cada cliente es asignado a uno nodo, y las restricciones (1.6i) son las restricciones de conservaci´on de flujo que tienen en cuenta que cada cliente absorbe una unidad de flujo del nodo 0 y otro del nodo n+1. Las restricciones Cap´ıtulo 1. Preliminares 13 (1.6j) imponen el requisito binario de que una arista e= (i, j)∈Eo es visitada por un anillo (suma de flujos igual a Q) o no (suma de flujos igual a 0). Las restricciones (1.6k) exigen la no negatividad de los flujos y las restricciones (1.6l) y (1.6m) exigen que las variables zij ywj tomen valor 0 o 1. En la Figura 1.3 podemos observar la distribuci´on del flujo que tendr´ıa uno de los anillos de la Figura 1.2. Para mayor claridad, se ha representado la r´eplica del nodo central ligeramente alejado de este. Se ejemplifica c´omo el flujo del camino que sale del nodo 0 es el n´umero de nodos que quedan por servir en ese anillo, y el que sale del nodo n+1 es la capacidad disponible del ciclo tras visitar cada nodo. Se puede comprobar que cada arista tiene flujo total Qunidades y que los nodos Steiner no ocupan capacidad en el anillo. Adem´as, se cumplen todas las restricciones de flujo. 0 2 5 8 9 11 12 13 16 17 19 22 25 27 11 5 10 6 9 7 7 9 6 10 5 11 4 12 2 14 016 Dep´osito Central Clientes Nodos Steiner Figura 1.3: Distribuci´on del flujo de uno de los anillos de la soluci´on del CmRSP de la Figura 1.2 con Q=16. 1.7 El problema multiobjetivo Como comentamos anteriormente, el problema RSP puede ser formulado como un problema biobjetivo en el que necesitamos minimizar la suma de los dos tipos de costes. As´ı, conservaremos las restricciones anteriores cambiando la funci´on objetivo por dos funciones por separado que debemos minimizar. Se vuelve a considerar un grafo dirigido G= (V, A) donde V es el conjunto de nodos y A el conjunto de arcos. Tambi´en usaremos las variables xij binarias para indicar si el arco (i,j) pertenece o no al ciclo y las variables yij binarias que son 1 si y solo si asignamos el nodo i al nodo j del ciclo. De nuevo, yii = 1 indica que el nodo i pertenece al ciclo. En general, podremos tomar la formulaci´on del RSP con aristas dirigidas. Pozo et al. (2024) presenta una metodolog´ıa para resolver problemas biobjetivo. En nuestro caso, tenemos los costes cij del ciclo y los costes dij de asignaci´on para cada arco 14 Secci´on 1.7. El problema multiobjetivo (i, j)∈A. Podemos definir la funci´on objetivo del problema como f(x) = (f1(x), f2(x)), donde f1(x) = X (i,j)∈A cijxij, f2(x) = X (i,j)∈A dijyij. Tomamos la definici´on usual de optimalidad o eficiencia de Pareto. Una soluci´on x∗ser´a ´optimo de Pareto si no existen soluciones ˆxque mejoren estrictamente alguno de los objetivos sin empeorar de forma simult´anea alguno de los otros, es decir, ∄ˆx:fi(ˆx)≤fi(x∗),∀i∈ {1,2}, y∃i∈ {1,2}:fi(ˆx)< fi(x∗). Si x∗es ´optimo de Pareto, f(x∗)∈R2ser´a un punto no dominado. Si fi(x∗)≤fi(ˆx),∀i∈ {1,2}, y ∃i∈ {1,2}:fi(x∗)< fi(ˆx), diremos que x∗domina a ˆxen el espacio de decisiones y f(x∗) domina a f(ˆx) en el espacio objetivo. As´ı, encontrar la soluci´on ´optima de un problema biobjetivo se reduce a encontrar el ´optimo de Pareto. En el siguiente cap´ıtulo se describe c´omo obtener los puntos no dominados (en el espacio objetivo), que pueden ser de dos tipos: soportados y no soportados. Los puntos no dominados soportados son aquellos que se pueden obtener al sumar dos valores objetivo, mientras que los no soportados no se pueden escribir de esta forma. Cap´ıtulo 2 M´etodos de resoluci´on En este cap´ıtulo se estudiar´an distintos m´etodos de resoluci´on para problemas de programaci´on lineal entera aplicables al bmRSP. En primer lugar, se desarrolla un m´etodo para la obtenci´on de los puntos no dominados, soportados y no soportados del problema. En segundo lugar, abordaremos el algoritmo de ramificaci´on, acotaci´on y corte que se emplea para simplificar la b´usqueda de soluci´on ´optima. Por ´ultimo, presentaremos un algoritmo heur´ıstico para el CmRSP que consta de dos fases: preasignaci´on y resoluci´on del CVRP. 2.1 M´etodo biobjetivo En el cap´ıtulo anterior se plante´o el bmRSP mediante una funci´on objetivo f(x) = (f1(x), f2(x)), donde f1(x) = X (i,j)∈A ci,jxi,j, f2(x) = X (i,j)∈A di,jyi,j. En Pozo et al. (2024) se proporciona un m´etodo para obtener el conjunto de puntos no dominados, soportados y no soportados; que se describe a continuaci´on. 2.1.1 Obtenci´on de los puntos no dominados soportados Con el objetivo de encontrar los puntos no dominados soportados, linealizamos la funci´on objetivo f(x) tomando λ∈[0,1] : 15 22 Secci´on 2.2. Algoritmo de ramificaci´on, acotaci´on y corte Z=maxkZk=max{18,21}= 21. Ahora bien, se tiene que el valor ´optimo es al menos 21 y la cota superior ¯ Z1= 20, por tanto, la soluci´on no puede estar en S1. As´ı, la rama S1del ´arbol de enumeraci´on sera podada por acotaci´on S S1S2 40 24 13 37 Max S S1S2 37 13 24 13 37 Figura 2.5: Ejemplo en el que podar no es posible Ejemplo 6. En la Figura 2.5, S vuelve a estar descompuesta en dos conjuntos S1yS2, indicando las cotas superior e inferior de cada uno. Se observa que ¯ Z=maxk¯ Zk= max{24,37}= 37 yZ=maxkZk=max{13,−∞} = 13. No se puede sacar ninguna conclusi´on y hay que seguir explorando ambos conjuntos. Bas´andonos en estos tres ejemplos, podemos dar al menos tres razones que nos permiten podar una rama y, de este modo, enumerar un gran n´umero de soluciones de forma impl´ıcita: 1. Podado por optimizaci´on: Zt={m´ax cx :x∈St}se ha encontrado y no es el m´aximo. 2. Podado por acotaci´on: ¯ Zt≤Z. 3. Podado por inviabilidad: St=∅. ¿C´omo conseguimos las cotas? La respuesta es sencilla, los l´ımites primarios (inferiores) los proporcionan soluciones factibles y los l´ımites duales (superiores) los obtendremos mediante relajaci´on o dualidad. La siguiente pregunta que debemos hacernos es c´omo construir un algoritmo de enumeraci´on impl´ıcito basado en las ideas anteriores. La Figura 2.6 nos da un diagrama de flujo con los pasos a seguir para encontrar la soluci´on ´optima, incluyendo un paso inicial heur´ıstico para comenzar desde un buen punto el algoritmo. Comenzamos sin cota inferior y sin una posible soluci´on x*. Tomando como inicializaci´on xHcuyo valor es ZH, tomaremos como primera cota inferior Z=ZHyx∗=xH. Cap´ıtulo 2. M´etodos de resoluci´on 23 Inicializaci´on: A˜nadir a una lista el problema inicial Scon politopo P. Tomamos Z=−∞ No se conoce x∗. Resolver con m´etodo heur´ıstico Tomar xHla soluci´on con valor ´optimo ZH Hacer Z=ZHyx∗=xH ¿La lista es vac´ıa? STOP Eliminar el problema Side la lista con la formulaci´on Piy la cota dual ¯ Zi Si ¯ Zi≤Z, podar por acotaci´on Resolver la relajaci´on lineal(LP) sobre Pi Tomar como cota dual ¯ Zi= valor objetivo de LP xi(LP) = soluci´on de LP Si Piest´a vac´ıo, podar por infactibilidad Si ¯ Zi≤Z, podar por acotaci´on Si xi(LP) es entero, actualizar la cota primal Z=¯ Ziyx∗=xi(LP). Podar por optimalidad Devolver dos subproblemas Si 1ySi 1 con politopos Pi 1yPi 2y cotas superiores ¯ Zi Si No No No No No S´ı S´ı S´ı S´ı S´ı Figura 2.6: Diagrama de flujo del algoritmo de ramificaci´on y acotaci´on. 24 Secci´on 2.2. Algoritmo de ramificaci´on, acotaci´on y corte 2.2.3 Planos de corte El cap´ıtulo 8 del libro de Wolsey (2020) nos da los algoritmos necesarios para incluir planos de corte en nuestro algoritmo B&B. Estos son los conocidos algoritmos de ramificaci´on, acotaci´on y corte (Branch &Cut -B&C). Primero, consideremos un problema general: (IP)max{cx :x∈X}(2.1) donde X={x:Ax ≤b, x ∈Zn +} Proposici´on 7. conv(X) = {x:˜ Ax ≤˜ b, x ≥0}es un poliedro. Este resultado nos permite reformular el problema 2.1 como un problema lineal: (LP)max{cx :˜ Ax ≤˜ b, x ≥0}(2.2) y as´ı, para cualquier valor de c, una soluci´on ´optima de punto extremo de (LP) es una soluci´on ´optima de (IP). El mismo resultado se cumple para problemas mixtos con X={(x, y)∈Zn +×Rp +:Ax +Gy ≤b}siempre que los datos A, G, b sean racionales. Desafortunadamente, para problemas NP-dif´ıciles casi no hay esperanza de encontrar una descripci´on “buena”de conv(x) . Dada una instancia de un problema NP-dif´ıcil, el objetivo en este apartado es encontrar formas efectivas de intentar aproximar su conv(X). Haremos uso del concepto de desigualdad v´alida. Recordemos que una desigualdad πx ≤π0es una desigualdad v´alida para X⊂Rnsi πx ≤π0∀x∈X. Sea X={x∈Zn:Ax ≤b}yconv(X) = {x:˜ Ax ≤˜ b}, las restricciones aix≤biy ˜aix≤˜ bison claramente desigualdades v´alidas para X. Debido a las numerosas desigualdades v´alidas que se pueden encontrar para X, nos plantearemos cu´ales de ellas son ´utiles para resolver una instancia concreta. Examinando la formulaci´on inicial P={x:Ax ≤b, c ≥0}con X=P∩Zn, encontraremos un conjunto de desigualdades v´alidas Qx ≤qpara X y las agregaremos a Pdando una nueva formulaci´on P′={x:Ax ≤b, Qx ≤q, x ≥0}con X=P′∩Zn. Aplicando el algoritmo B&Ca la formulaci´on P′, surgen varias ventajas y desventajas:  Ventajas: El algoritmo se vuelve m´as eficiente. Si las desigualdades v´alidas est´an bien elegidas, de modo que la formulaci´on P’ sea significativamente m´as peque˜na que P, las cotas deber´ıan mejorar y, por lo tanto, el algoritmo de B&Cdeber´ıa ser m´as efectivo. Adem´as, aumentan las posibilidades de encontrar soluciones enteras factibles en el transcurso del algoritmo. Cap´ıtulo 2. M´etodos de resoluci´on 25  Desventajas: A menudo, la familia de desigualdades v´alidas que nos gustar´ıa agregar es demasiado grande. En tales casos, o bien la formulaci´on se vuelve muy extensa y requiere de mucho tiempo encontrar la soluci´on, o se vuelve imposible utilizar el algoritmo de B&Cdebido a las numerosas restricciones. Para buscar desigualdades b´asicas, observamos que, en muchas ocasiones, la regi´on factible Xpuede ser escrita naturalmente como la intersecci´on de dos o m´as conjuntos, es decir, X=X1∩X2. Descomponer o centrarse en uno de los dos conjuntos puede ser una buena idea. Si el problema de optimizaci´on sobre el conjunto X2=P2∩Znes f´acil de resolver, se puede encontrar una descripci´on expl´ıcita de P′2=conv(P2∩Zn). En este caso, conviene reemplazar la formulaci´on inicial P1∩P2por una formulaci´on mejorada P′=P1∩P′2. Incluso si el problema sobre X2no fuera f´acil de resolver, podr´ıamos encontrar desigualdades v´alidas para X2que nos ayudar´ıan a mejorar su formulaci´on a P′2⊂P2. As´ı, se proporciona nuevamente una formulaci´on mejorada P′=P1∩P′2para X. Algoritmo 4: Un algoritmo de planos de corte Entrada: Establecer t= 0 y P0=P; 1while no se cumpla la condici´on de parada do 2Iteraci´on t: Resolver el programa lineal: ¯ Zt= m´ax{cx |x∈Pt} Sea xtuna soluci´on ´optima; 3if xt∈Znthen 4STOP; 5xtes una soluci´on ´optima para el problema. 6else 7Resolver el problema de separaci´on para xty la familia F. 8if encontramos una desigualdad (πt, πt 0)∈ F con πtx>πt 0que corta a xtthen 9Tomamos Pt+1 =Pt∩ {x:πtx≤πt 0}y aumentamos t. 10 else 11 STOP Ahora bien, consideremos X=P∩ZnyFuna familia de desigualdad v´alidas donde (π, π0)∈ F para Xsi πx ≤π0. En numerosos casos, Fcontiene demasiadas restricciones (2n o m´as) como para ser a˜nadidas a priori. Adem´as, dada una funci´on objetivo en espec´ıfico, el objetivo no es encontrar la envolvente convexa completa sino encontrar una buena aproximaci´on 26 Secci´on 2.3. Algoritmo heur´ıstico de dos fases de conv(X) en el plano de una soluci´on ´optima. El Algoritmo 4 describe un algoritmo b´asico de planos de corte para (IP)max{cx :x∈ X}que genera desigualdades ´utiles a partir de F. Si el algoritmo termina sin haber encontrado una soluci´on entera para (IP), Pt=P∩ {x:πix≤πi 0, i = 1, . . . , t} es una formulaci´on mejorada que puede ser la entrada de un algoritmo B&C. Cabe se˜nalar que en la pr´actica es mejor a˜nadir varios cortes en cada iteraci´on en vez de uno a la vez. 2.3 Algoritmo heur´ıstico de dos fases Baldacci et al. (2007) presenta un algoritmo heur´ıstico de 2 fases para resolver el CmRSP. La primera consiste en obtener un preassignment que determine la asignaci´on de cada nodo. La segunda fase toma los nodos que se han asignado a s´ı mismos y resuelve el CVRP (ver Cap´ıtulo 1) formado por ese conjunto. En el art´ıculo, esa fase es llamada heur´ıstica primal. Algoritmo 5: Algoritmo de preasignaci´on(y∗, σ, S, q) Entrada: Variables y∗. Salida: Asignaciones σ, conjunto S de nodos con clientes asignados y sus pesos q. 1for i∈Udo 2σ(i) := −1; 3¯ U:= U,clientes que tienen que ser asignados ; 4N:= V′,nodos a los que un cliente puede ser asignado ; 5S:= ∅,nodos usados para asignar; 6Preasignaci´on de clientes 7while ¯ U=∅do 8Selecciona aleatoriamente a un cliente i∈¯ U; 9¯ U:= ¯ U\ {i}; 10 σ(i) = argmax{y∗ ij :j∈Ci∩N}; 11 S:= S∪ {σ(i)}; 12 if σ(i)=ithen 13 N:= N\ {i} 14 for j∈Sdo 15 qj:= |{i∈U:σ(i) = j}| En primer lugar, para obtener el preassigment, resolvemos la relajaci´on lineal del problema Cap´ıtulo 2. M´etodos de resoluci´on 27 obteniendo las variables y* de la soluci´on fraccional. Asignaremos a cada nodo i, el nodo j que tenga mayor y∗ ij. El procedimiento de preasignaci´on (Algoritmo 5) nos devuelve: 1. El nodo σ(i) que asignaremos a cada nodo i∈U. 2. El conjunto S⊂V′de los nodos que tienen algunos clientes asignados. 3. El n´umero q(j) de clientes asignados a cada nodo j∈S. En segundo lugar, el algoritmo heur´ıstico resuelve el CVRP definido por el conjunto de clientes en S con las demandas q y con m veh´ıculos de capacidad Q. Con esta soluci´on se construir´a la soluci´on del CmRSP a˜nadiendo a la del CVRP las asignaciones de la fase de preasignaci´on. 28 Secci´on 2.3. Algoritmo heur´ıstico de dos fases Cap´ıtulo 3 Implemetaci´on Una vez estudiada la modelizaci´on del problema y los m´etodos de resoluci´on que pueden ayudarnos a encontrar la soluci´on, el siguiente paso ser´a aplicar estas t´ecnicas para encontrar la soluci´on de un conjunto de instancias del problema. Para ello, se utilizar´an programas inform´aticos donde se implementar´an los diferentes modelos, algoritmos y elementos de an´alisis de las soluciones obtenidas. En este cap´ıtulo se abordar´a la implementaci´on en Python de los modelos matem´aticos vistos en el Cap´ıtulo 1 para formular el bmRSP. Tambi´en se implementar´an los algoritmos de resoluci´on del Cap´ıtulo 2. En primer lugar, se explicar´a el c´odigo usado para generar instancias que nos servir´an como datos de entrada en modelos y algoritmos. Posteriormente, se analizar´a la implementaci´on de los 3 algoritmos que se vieron en el cap´ıtulo anterior: el algoritmo de ramificaci´on, acotaci´on y corte; el algoritmo heur´ıstico de 2 fases y el algoritmo de la ε-restricci´on. Se mostrar´an los resultados obtenidos con cada uno de ellos, junto con representaciones gr´aficas y tablas de comparaci´on entre instancias. 3.1 Generaci´on de instancias Una instancia del bmRSP es un conjunto de datos para el problema que ser´a guardado en un archivo de texto. Se han tomado grafos completos siendo V={1, ..., n}el conjunto de nodos con coordenadas enteras del 1 al 10. Tanto los costes del anillo cij como los de asignaci´on dij ser´an las distancias entre los nodos extremos del arco (i, j). Para crear distintas instancias con este formato, implementaremos un generador de instancias en Python que redacte un archivo de texto distinto seg´un el tama˜no y el n´umero de anillos que se le proporcione. Indicar´a en la primera l´ınea el n´umero de anillos mdebe tener la soluci´on del bmRSP. Tambi´en se incluir´a 29 30 Secci´on 3.2. Carga de datos y representaci´on de soluciones el n´umero de nodos ndel grafo y generar´a sus coordenadas usando la librer´ıa random dentro del cuadrado (0,10) ×(0,10). El generador de instancias duplicar´a el nodo 0 con el fin de contar con la r´eplica del dep´osito central en el nodo n+ 1. A˜nadir´a las aristas y arcos del grafo completo de n+ 1 nodos (ya que contaremos tambi´en con implementaciones no dirigidas que har´an uso de las aristas) y los cardinales de estos conjuntos que ser´an necesarios a la hora de cargar las instancias. Generaremos tambi´en la lista de arcos y aristas junto a sus costes que ser´an determinados por la distancia eucl´ıdea que hay entre sus extremos. Por simplicidad y con el fin de acotar todas las posibilidades que ofrecen los experimentos, se ha optado por definir el conjunto de nodos Steiner como vac´ıo (W=∅), por tanto, V′=UyV=U∪ {0, n + 1}son los ´unicos conjuntos con los que trabajaremos en nuestras implementaciones. Se generan 3 instancias de cada tama˜no con el objetivo de comparar resultados y obtener los valores medios de los par´ametros que se obtienen en la resoluci´on de los algoritmos (valor objetivo, tiempo de resoluci´on, gaps, etc). 3.2 Carga de datos y representaci´on de soluciones En todos los modelos y algoritmos que veremos a continuaci´on, el primer paso ser´a hacer una lectura de los datos del archivo de texto que hemos generado. Esa lectura nos devolver´a el n´umero de nodos, el n´umero de anillos y los distintos costes del problema. Para ello, se crea una funci´on de lectura que abre el archivo de texto que le demos como argumento de entrada, lee todas las l´ıneas del archivo y las guarda en una lista. Tomar´a los valores de m,n, el n´umero de arcos y el de aristas del grafo. Se inicializan listas vac´ıas “aristas” y “arcos” para almacenarlas. Se recorren las l´ıneas correspondientes a aristas y arcos, se eliminan los caracteres de separaci´on, se convierten en tuplas de n´umeros enteros utilizando las funciones map, tuple, strip ysplit y se almacenan en las listas anteriormente mencionadas. De manera parecida, almacenamos en diccionarios “costes ciclos” y “costes asignaciones” las l´ıneas correspondientes de la instancia. Si se pretende realizar una representaci´on gr´afica de la soluci´on que vayamos a obtener, la funci´on tambi´en genera un c´odigo LaTeX llamado “contenido latex” que incluye el encabezado del archivo e inicia un entorno tikzpicture para dibujar el gr´afico. Luego, se utiliza un bucle for para iterar sobre las l´ıneas que representan los nodos en el archivo de instancia. Por cada nodo, se extraen sus coordenadas y se genera un comando LaTeX para dibujar el nodo en esas coordenadas con un n´umero asociado. Tambi´en se a˜naden los comandos que dibujan los resultados del algoritmo correspondiente. Para ello, se crean bucles for que recorren los valores de las distintas variables, tomando las que est´an activadas y a˜nadiendo para cada una de ellas la correspondiente l´ınea de comando que dibuja ese arco. Normalmente, los arcos se dibujaran Cap´ıtulo 3. Implemetaci´on 31 como flechas ya que la soluci´on es dirigida. Sin embargo, en el caso del problema de flujo los valores positivos de fse dibujar´an como aristas ya que este problema considera los ciclos no dirigidos. Se ha a˜nadido tambi´en una l´ınea de comando que muestra en el dibujo el coste total de la soluci´on obtenida. Se a˜nadir´a el cierre del entorno tikzpicture y del documento y, por ´ultimo, se generar´a un archivo de LaTeX con el contenido completo. 3.3 Descripci´on general de la resoluci´on Aplicando cualquier m´etodo de resoluci´on, encontrar la soluci´on del problema puede ser computacionalmente muy costoso. Como este proceso podr´ıa alargarse demasiado en el tiempo, le daremos un time-limit al algoritmo. Los elementos reportados en las tablas de este cap´ıtulo son: 1. ObjUB: La mejor soluci´on factible que haya encontrado el solver hasta el momento. Si el tiempo de resoluci´on es menor que el tiempo l´ımite que le hemos indicado, ser´a la soluci´on ´optima. 2. ObjLR: El valor de la relajaci´on lineal del problema antes de empezar a ramificar. 3. ObjLB: El valor de la cota inferior que tiene el algoritmo cuando ha llegado al tiempo l´ımite. Si el tiempo de resoluci´on es menor que el tiempo l´ımite que le hemos indicado, ser´a la soluci´on ´optima. 4. time: El tiempo de ejecuci´on que ha invertido el modelo en encontrar la soluci´on ´optima factible, que ser´a igual al time-limit en caso de no haber podido encontrarla. 5. nod: El n´umero de nodos explorados en el ´arbol de busqueda al aplicar el algoritmo B&B. 6. gUB: El GAP de optimalidad gUB =ObjUB −ObjLB ObjUB ∗100 %. 7. gUR: El GAP de la relajaci´on lineal gUR =ObjUB −ObjLR ObjUB ∗100 %. En la siguiente secci´on, se resuelve el CmRSP (es decir, el bmRSP con un ´unico objetivo). Posteriormente, se resolver´a el bmRSP en el apartado 3.5. 3.4 Algoritmo de ramificaci´on, acotaci´on y corte Python-Gurobi es una potente interfaz que permite a los usuarios aprovechar las capacidades de optimizaci´on de Gurobi, uno de los solver l´ıderes en el ´ambito de la programaci´on 38 Secci´on 3.5. Algoritmo heur´ıstico de dos fases 1. El modelo de flujo no admite formulaci´on en grafos dirigidos. En t´erminos de valor objetivo, esto no es relevante si los costes son determinados por las distancias entre los nodos. Sin embargo, si se toman costes no sim´etricos en los que xij =xji el modelo de flujo no ser´ıa v´alido. Se muestra en la Figura 3.3 la soluci´on de una instancia de 50 nodos que el modelo de eliminaci´on de subtours obtuvo a los 284 segundos y el modelo de flujo s´olo necesit´o 14 segundos. Se puede observar una peque˜na diferencia en el nodo 44 debido a que no cambia el coste final, por lo que ambas son soluciones ´optimas y podr´ıan haber salido en cualquiera de los modelos. 2. El modelo de eliminaci´on de subtours usa m´as variables binarias porque admite grafos dirigidos. 3. El modelo de flujo no se puede modificar para admitir soluciones con menos de manillos. Para saber cu´al es la soluci´on ´optima con menos de manillos tendr´ıamos que resolver el problema para m′= 1, . . . , m y comparar los valores objetivo. 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1617 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50510 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1617 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 5051 Funci´on objetivo: 67.21 (a) Producida por el modelo de eliminaci´on de subtours y el algoritmo B&C 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1617 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50510 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1617 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 5051 Funci´on objetivo: 67.21 (b) Producida por el modelo de flujo Figura 3.3: Soluciones de la 2 ª instancia de 50 nodos y Q=33 obtenidas con distintos modelos. 3.5 Algoritmo heur´ıstico de dos fases Como vimos en la Secci´on 2, existe un algoritmo heur´ıstico para el CmRSP que consta de 2 fases: preasignaci´on y resoluci´on del CVRP. Recordemos que el algoritmo de preasignaci´on Cap´ıtulo 3. Implemetaci´on 39 se basaba en la soluci´on de la relajaci´on lineal del problema. Por tanto, para la implementaci´on del algoritmo en Python debemos calcular unos valores iniciales y∗. Para ello, redefinimos xey como variables continuas dentro del dominio x, y ∈[0,1]. Llamaremos a la funci´on que resuelve el modelo relajado y guardaremos las variables de asignaci´on de la soluci´on obtenida. Como nuestras instancias no incluyen nodos Steiner, trabajaremos con U=V′yV=U∪ {0, n + 1}. Se inicializan varias estructuras de datos:  ¯ U:= U, conjunto de clientes que tienen que ser asignados (en la inicializaci´on faltan todos por asignar).  N:= U, conjunto de nodos a los que un cliente puede ser asignado.  S:= [ ], lista que se usar´a para guardar los nodos a los que se les asignan otros.  σ=i:−1/i ∈U, inicializamos con este diccionario que le asigna a cada nodo el −1. Creamos un bucle que mientras ¯ Utenga elementos, coja alguno de ellos elimin´andolo del conjunto, le asigne el nodo de Nque maximice el valor de la variable de asignaci´on que les corresponde en la soluci´on y∗del problema relajado. As´ı, el nodo elegido se a˜nade al conjunto S. Si el cliente es asignado a un nodo diferente del suyo, se elimina de la lista N de nodos disponibles. Cuando terminamos de asignar todos los nodos de ¯ U,σrecoge todas las preasignaciones y Stodos los nodos que deben pertenecer a alg´un anillo. Para terminar la primera fase, se calculan los pesos qde cada nodo basados en el n´umero de clientes asignados a cada uno. Como este algoritmo no incluye al dep´osito central, lo a˜nadimos posteriormente a Sya que aunque no se le puede asignar nodos, s´ı debe estar en los anillos de la segunda fase. Como para la siguiente fase necesitamos guardar las distancias de la instancia y los nodos, aprovechamos la lectura de los datos para crear “contenido latex” que guarda las lineas de comando necesarias para dibujar la soluci´on de la preasignaci´on. A modo de ejemplo, una instancia de 25 nodos devolver´ıa tras la primera fase las preasignaciones que se muestran en la Figura 3.4a. La segunda fase consiste en resolver el problema CVRP definido por los nodos del conjunto S con mveh´ıculos, con el mismo dep´osito central, la capacidad Qy tomando q(i) (n´umero de asignaciones que llegan al nodo i) como la demanda del cliente i. Por la estructura del algoritmo de preasignaci´on, Spuede contener repeticiones y las eliminamos convirti´endolo en un conjunto. Cargamos la instancia para tomar los costes de ciclo que hay entre los nodos de S. Se crea el modelo de rutas que nos proporciona Leite (2023) utilizando la biblioteca Google OR- Tools que usa m´etodos de b´usqueda local asociados a algoritmos constructivos. Para simplificar la implementaci´on, creamos un diccionario que mapea cada nodo de S con su enumeraci´on, as´ı conseguimos transformar S→1, . . . , n2donde n2=|S|. Se configuran los callbacks para calcular las distancias y las demandas creando nuevos diccionarios con el renombramiento de 40 Secci´on 3.6. Algoritmo ´epsilon-restricci´on los clientes. Una vez tenemos todo lo necesario, se establecen los par´ametros de b´usqueda del solver, incluyendo la estrategia de la primera soluci´on y la metaheur´ıstica de b´usqueda local. Para ello, usamos:  First solution strategy: Esta configuraci´on determina la estrategia que se utilizar´a para encontrar la primera soluci´on del problema de rutas. Seleccionamos CHRISTOFIDES que se refiere al algoritmo de Christofides,conocido por ser eficiente para encontrar soluciones iniciales para problemas de rutas.  Local search metaheuristic: Esta configuraci´on determina la metaheur´ıstica de b´usqueda local que se utilizar´a para mejorar las soluciones encontradas. Seleccionamos GUIDED LOCAL SEARCH, que es una de las metaheur´ısticas disponibles en OR-Tools. La b´usqueda local guiada es un enfoque heur´ıstico que explora el espacio de soluciones al realizar peque˜nos cambios en una soluci´on actual y evaluando si estos cambios conducen a mejoras.  Time limit: Esta configuraci´on establece un tiempo l´ımite para la b´usqueda de soluciones. En este caso, se ha establecido en 300 segundos (5 minutos). Esto significa que el algoritmo de b´usqueda intentar´a encontrar una soluci´on dentro de este l´ımite de tiempo. Si no se encuentra, la b´usqueda se detendr´a y se devolver´a la mejor soluci´on encontrada hasta ese momento. Por ´ultimo, se resuelve el problema del CVRP utilizando el modelo de rutas configurado y los par´ametros de b´usqueda especificados. Finalizamos almacenando en listas el recorrido de cada veh´ıculo y accediendo al coste de ciclos total. La soluci´on obtenida aplicando este algoritmo al conjunto de nodos Ses la Figura 3.4b. Resulta evidente que la soluci´on encontrada no es la ´optima, sin embargo, esto es natural al tratarse de un algoritmo metaheur´ıstico. En la Tabla 3.5, podemos observar las soluciones que proporciona el algoritmo de dos fases. Se han tomado las cotas inferiores y el valor de la relajaci´on lineal de la tabla de flujo. Se observa que aunque el algoritmo es muy r´apido, no llega a soluciones ´optimas. Usando la soluci´on del CVRP y las asignaciones de la primera fase, obtenemos la soluci´on del CmRSP que podemos ver en la Figura 3.4c. Al tomar como demanda de cada nodo el n´umero de preasignaciones que se le asignaron en la primera fase, la construcci´on del CVRP hace que se cumpla la capacidad de cada anillo en la soluci´on final. 3.6 Algoritmo ´epsilon-restricci´on Para implementar en Python el algoritmo ´epsilon-restricci´on que proporciona el conjunto de soluciones no dominadas del bmRSP, necesitamos modelar el problema con f2como funci´on Cap´ıtulo 3. Implemetaci´on 41 |V|Q inst gUR gUB time ObjLR ObjLB ObjUB 10 6 1 20.71 10.05 0.07 30.50 34.60 38.47 10 6 2 59.45 42.12 0.05 29.76 42.52 73.45 10 6 3 53.28 44.70 0.05 42.13 49.85 90.15 50 33 1 61.88 56.89 0.82 60.44 68.31 158.52 50 33 2 62.54 55.56 0.75 56.68 67.21 151.25 50 33 3 50.13 43.68 0.66 62.34 70.40 125.01 100 66 1 57.50 56.23 7.52 102.79 105.88 241.89 100 66 2 62.14 61.83 6.94 102.00 102.82 269.45 100 66 3 62.83 62.61 7.57 102.18 102.82 274.93 200 133 1 75.30 73.24 118.08 217.69 235.84 881.77 200 133 2 74.70 73.51 122.89 216.28 226.66 855.37 200 133 3 75.67 74.15 118.11 221.77 235.58 911.31 Tabla 3.5: Soluciones obtenidas con el algoritmo heur´ıstico de dos fases para m=2. 0 1 2 3 4 5 6 7 8 9 10 11 Nodo: peso 3:2, 5:2, 2:1, 1:3, 10:1, 7:1 0 1 2 3 5 7 10 0 (a) Primera fase y pesos del anillo 0 1 2 3 4 5 6 7 8 9 10 11 (b) Segunda fase 0 1 2 3 4 5 6 7 8 9 10 11 (c) Soluci´on final Figura 3.4: Resoluci´on por pasos del algoritmo heur´ıstico en una instancia de 10 nodos, 2 anillos y Q=10 42 Secci´on 3.6. Algoritmo ´epsilon-restricci´on objetivo. Es decir, creamos un modelo con Gurobi a˜nadiendo todas las restricciones del CmRSP y minimizando ´unicamente los costes de asignaci´on. Ese modelo nos devolver´a un valor ´optimo Z1 2y podremos calcular los costes de anillo Z1 1con los valores que se le han dado a las variables. El punto no dominado (Z1 1, Z1 2) ser´a la inicializaci´on de nuestro algoritmo. Crearemos una lista [(Z1 1, Z1 2)] que a˜nade el punto no dominado que nos devuelve cada paso del algoritmo. La idea es imponer que el coste de los anillos decrezca en cada iteraci´on con respecto a la anterior al menos un valor ε. Esto lo conseguimos a˜nadiendo en la iteraci´on ila restricci´on: X (i,j)∈A cijxij ≤Zi−1 1−ε(3.1) Se crea un bucle for para realizar las iteraciones i∈ {1, ..., 100}, asegurando que el algoritmo pare en un tiempo considerable aunque no haya llegado a minimizar en su totalidad los costes del ciclo. Al minimizar el coste de asignaci´on a˜nadiendo la restricci´on (3.1), estaremos obteniendo un nuevo punto no dominado (Zi 1, Zi 2). Llegado el momento en el que no sea posible minimizar los costes del ciclo, el modelo se volver´a infactible y parar´a el bucle. Por ello, en nuestro algoritmo a˜nadimos un bloque try-except que comprueba la factibilidad y, en caso de ser necesario, imprime un mensaje diciendo en qu´e iteraci´on se ha vuelto infactible y termina el proceso. Por ´ultimo, imprimimos la lista de puntos no dominados que hemos encontrado en cada iteraci´on. Podemos observar c´omo evoluciona la soluci´on del problema en cada iteraci´on con una instancia de 10 nodos y dos anillos. El primer paso, minimiza ´unicamente el coste de asignaci´on, por tanto, es natural que los anillos formados no sean ´optimos. Este fen´omeno es visible en la Figura 3.5. 0 1 2 3 4 5 6 7 8 9 10 110 1 2 3 4 5 6 7 8 9 10 11 Coste de ciclo: 60.83, Coste de asignaci´on: 0.0 Figura 3.5: Primera paso: minimizar el coste de asignaci´on. Ahora bien, como en cada iteraci´on disminuye el coste de ciclo, se ir´an consiguiendo mejoras en los anillos y a˜nadiendo asignaciones a la soluci´on. Esto podemos verlo en las Figuras Cap´ıtulo 3. Implemetaci´on 43 3.6. 0 1 2 3 4 5 6 7 8 9 10 110 1 2 3 4 5 6 7 8 9 10 11 Coste de ciclo: 38.62, Coste de asignaci´on: 0.0 (a) Iteraci´on 4 0 1 2 3 4 5 6 7 8 9 10 110 1 2 3 4 5 6 7 8 9 10 11 Coste de ciclo: 32.99, Coste de asignaci´on: 2.24 (b) Iteraci´on 10 0 1 2 3 4 5 6 7 8 9 10 110 1 2 3 4 5 6 7 8 9 10 11 Coste de ciclo: 27.19, Coste de asignaci´on: 7.4 (c) Iteraci´on 20 0 1 2 3 4 5 6 7 8 9 10 110 1 2 3 4 5 6 7 8 9 10 11 Coste de ciclo: 22.11, Coste de asignaci´on: 12.77 (d) Iteraci´on 30 Figura 3.6: Soluciones obtenidas en diferentes iteraciones del algoritmo de la ε-restricci´on. En la ´ultima iteraci´on, el coste de ciclo no puede minimizarse m´as por lo que devolver´a una soluci´on con dos anillos m´ınimos y el resto de nodos se incluir´an a la red a trav´es de asignaciones. Esto podemos verlo en las Figuras 3.7. 0 1 2 3 4 5 6 7 8 9 10 110 1 2 3 4 5 6 7 8 9 10 11 Coste de ciclo: 6.48, Coste de asignaci´on: 31.28 Figura 3.7: ´ Ultima iteraci´on. Este m´etodo puede resultar ´util si tenemos mayor inter´es en minimizar uno de los dos 44 Secci´on 3.6. Algoritmo ´epsilon-restricci´on costes, pudiendo encontrar el balance que m´as se ajuste a nuestras necesidades. Debemos tener en cuenta que la funci´on objetivo del problema siempre es el coste de asignaci´on, por lo que no nos asegura que los anillos sean ´optimos. En la Figura 3.6a se observa que ha mejorado el coste del ciclo respecto de la Figura 3.5 y, sin embargo, el ciclo (0,6,1,10,3,11) ser´ıa de menor longitud si tom´aramos la reordenaci´on (0,3,6,10,1,11). 5 10 15 20 25 30 35 0 10 20 30 9 : 34,91 11 13 14 : 34,28 15 16 17 19 : 34,38 20 22 : 34,59 23 : 34,7 24 : 35,33 25 26 27 : 34,98 28 29 : 35,61 30 : 34,81 31 32 : 34,88 33 : 35,2 34 : 34,98 37 38 : 36,52 39 : 35,2 40 41 : 36,67 43 : 36,76 45 : 37,76 46 48 49 : 37,76 z1 z2 Puntos no dominados Figura 3.8: Lista de puntos no dominados que nos devuelve el algoritmo de la ε-restricci´on en una instancia de 10 nodos. Si quisi´eramos encontrar el coste total m´ınimo, tendr´ıamos que buscar en la lista de pares que nos devuelve el algoritmo, aquel con suma m´ınima y seleccionar la soluci´on de la iteraci´on a la que corresponde. En este caso, hemos representado en la Figura 3.8 cada punto no dominado con una etiqueta de la forma Iteraci´on: Suma de costes. Se han omitido algunas sumas irrelevantes para mayor claridad. N´otese que los valores objetivo de algunas iteraciones no est´an representados. Esto es debido a que son puntos dominados, los cuales elimina el algoritmo en el ´ultimo paso. Se observa que la iteraci´on que genera un grafo de menor coste es la Iteraci´on 14, con un coste de 34,28 unidades. Cap´ıtulo 4 Ap´endice 4.1 Formulaciones no dirigidas del RSP y el CmRSP Labb´e et al. (2004) formula el RSP en grafos no dirigidos como un programa lineal entero mixto. Plantea el problema con un grafo G= (V, A ∪E) donde V={1, . . . , n}es el conjunto de v´ertices (consideraremos n≥5 para evitar casos triviales), E={[i, j] : i, j ∈V, i < j}es el conjunto de aristas y A={(i, j) : i, j ∈V}es el conjunto de arcos incluyendo los bucles (i, i). Podemos observar que ahora los arcos son dirigidos pero las aristas no, por tanto,tendremos la mitad de variables relacionadas con el anillo (en consecuencia, habr´a menos restricciones). Cada arista [i, j]∈Etendr´a un coste de anillo cij y cada arco (i, j)∈Atendr´a un coste de asignaci´on dij. Por tanto, para solucionar el RSP hay que encontrar un ciclo en un subconjunto de v´ertices V′∈Vque incluya 1 (que se denotar´a ra´ız o dep´osito) y, al menos, dos v´ertices m´as; minimizando el coste total. A cada v´ertice ino incluido en el anillo asignaremos uno jde forma que cumpla que dij =minvk∈V′dik. As´ı, el coste que habr´a que minimizar ser´a la suma de los costes del ciclo y los de asignaci´on que ser´an Pi∈V\V′minj∈V′dij. Tomaremos variables binarias xij que toman el valor 1 si la arista [i, j]∈Eforma parte del ciclo y 0 en caso contrario. Por otro lado, para cada arco (i, j)∈Ase toma la variable yij que vale 1 si ies asignado a un v´ertice jdel ciclo. Ahora, para modelarlo necesitamos definir los conjuntos E(S) := {[i, j]∈E:i, j ∈S}yδ(S) := {[i, j],[j, i]∈E:i∈S, j /∈S}para cada S⊂V. Denotaremos δ(i) a las aristas de icon el resto de v´ertices jcon i<j. Para cada E′⊂Edefiniremos x(E′) := P[i,j]∈E′xij. As´ı, el problema ser´a 45 46 Secci´on 4.1. Formulaciones no dirigidas del RSP y el CmRSP m´ın X [i,j]∈E cijxij +X (i,j)∈A dijyij (4.1a) s.a. x(δ(i)) = 2yii ∀i∈V, (4.1b) X j∈V yij = 1 ∀i∈V\ {1},(4.1c) x(δ(S)) ≥2X j∈S yij = 1 ∀S⊂V: 1 /∈S, i ∈S, (4.1d) xij ∈ {0,1} ∀[i, j]∈E, (4.1e) yij ≥0∀(i, j)∈E, (4.1f) y11 = 1 (4.1g) y1j= 0 ∀j∈V\ {1},(4.1h) yjjentero ∀j∈V\ {1},(4.1i) De esta forma, conseguimos minimizar la suma de los distintos costes imponiendo las condiciones necesarias para tener una soluci´on factible. Como x(δ(i)) denota el n´umero de aristas del ciclo que tienen a icomo v´ertice y la variable yii ser´a 1 si y solo si ipertenece al anillo, la condici´on (4.1b) asegura que el grado de un v´ertice sea 2 si y solo si pertenece al ciclo. Por otro lado, la condici´on (4.1c) asegura que todos los v´ertices que no sean el dep´osito (que tendr´a que ser parte del ciclo obligatoriamente) tendr´an 2 opciones: formar parte del ciclo o ser asignados a un v´ertice del ciclo. La restricci´on (4.1d) impone la conectividad, obligando a que todo conjunto S debe estar conectado a su complementario al menos por dos aristas del ciclo a no ser que uno de los v´ertices vi∈Sest´e asignado a vj∈S. La combinaci´on de (4.1b),(4.1e),(4.1g) y (4.1h) garantiza que la soluci´on incluya al menos un ciclo conteniendo al dep´osito. Como queremos que haya un ´unico ciclo, lo imponemos en la restricci´on (4.1d). Tambi´en imponemos que los v´ertices que no pertenezcan al ciclo se asignen a uno que s´ı con las restricciones (4.1b),(4.1c),(4.1d) y (4.1f). Imponemos la condici´on de integrabilidad ´unicamente para las yjj porque para las variables yij(vi=vj) no es necesario, ya que la funci´on objetivo se minimizar´a cuando asignemos los v´ertices que no est´an en el ciclo al m´as cercano del anillo. Impl´ıcitamente, hemos garantizado que el ciclo incluya al menos tres v´ertices incluyendo la ra´ız y que la comunicaci´on entre v´ertices permanezca a´un fallando alguno de los del anillo. La formulaci´on del CmRSP en un grafo G= (V, E ∪A) no dirigido sigue la misma notaci´on que en la Secci´on 1.6 a˜nadiendo el conjunto Ede aristas del grafo que esta vez definiremos como E={(i, j) : i, j ∈V, i =j}. Para cada E′⊂E, denotaremos como V(E′) al conjunto de v´ertices incidentes con al menos una arista de E′. Recordemos que V=V′∪0, n + 1 Cap´ıtulo 4. Ap´endice 47 son los nodos del grafo incluyendo el dep´osito central y su copia y V′=U∪Wdonde Ues el conjunto de clientes y Wel conjunto de nodos Steiner. Las conexiones entre dos v´ertices(nodos incluidos en alg´un anillo) ser´an aristas y las conexiones directas entre nodos y v´ertices ser´an arcos. Se extender´an algunas sumas a los pares (i, j)/∈Aasumiendo que entonces las sumas adicionales valen 0. Como hicimos en la formulaci´on del RSP, denotaremos como δ(S) = {(i, j)∈E:i∈ S, j /∈S}al conjunto de aristas que tengan un ´unico nodo en S, para todo S⊂V; y δ(i) en el caso de que S={i}. Tambi´en denotaremos U(S) = S∩Ual conjunto de clientes en S. m´ın X (vi,vj)∈E cijxij +X (vi,vj)∈A dijyij (4.2a) s.a. X (i,j)∈δ(0) xij =m, (4.2b) X (i,j)∈δ(n+1) xij =m, (4.2c) X (i,j)∈δ(i) xij = 2yii ∀i∈U, (4.2d) X (i,j)∈δ(i) xij = 2wi∀i∈W, (4.2e) X j∈V′ yij = 1 ∀i∈U, (4.2f) X (i,j)∈δ(S) xij ≥2 QX i∈UX j∈S yij ∀S⊆V′:S=∅(4.2g) xij ∈ {0,1} ∀(i, j)∈E, (4.2h) yij ∈ {0,1} ∀(i, j)∈A, (4.2i) wi∈ {0,1} ∀i∈W, (4.2j) La expresi´on (4.2a) representa la funci´on objetivo que trata de minimizar el coste total obtenido como la suma de los costes de conexi´on de los arcos de los anillos y los costes de conexi´on directa de los clientes que se conectan a los nodos de los anillos. Las restricciones (4.2b) y (4.2c) aseguran que m anillos salen del nodo 0 y entran por el nodo n + 1. Las restricciones (4.2d) y (4.2e) imponen que el grado de cada nodo i∈V′es 2 si pertenece a un anillo. Estas dos restricciones permitir´ıan eliminar una de las restricciones (4.2b) o (4.2c) por ser redundantes, pero se mantienen para simplifcar la presentaci´on. Las restricciones (4.2f) establecen que un cliente i o bien est´a en un anillo (yii = 1) o est´a asignado a un nodo de alg´un anillo (yij = 1 para alg´un j∈Ci). Las restricciones (4.2g) son las restricciones fraccionales de 54 Bibliograf´ıa Baldacci, R. and Dell’Amico, M. (2010) “Heuristic algorithms for the multi-depot ring-star problem”. European Journal of Operational Research, vol. 203(1): 270–281. Baldacci, R., Dell’Amico, M. and Gonz´alez, J. S. (2007) “The capacitated m-ring-star problem”. Operations research, vol. 55(6): 1147–1162. Baldacci, R., Hill, A., Hoshino, E. A. and Lim, A. (2017) “Pricing strategies for capacitated ring-star problems based on dynamic programming algorithms”. European Journal of Operational Research, vol. 262(3): 879–893. Calvete, H. I., Gal´e, C. and Iranzo, J. A. (2013a) “An efficient evolutionary algorithm for the ring star problem”. European Journal of Operational Research, vol. 231(1): 22–33. Calvete, H. I., Gal´e, C. and Iranzo, J. A. (2013b) “An evolutionary algorithm for the biobjective capacitated m-ring star problem”. In P. Perny, M. Pirlot and A. Tsouki`as, eds., Algorithmic Decision Theory, pp. 116–129. Springer Berlin Heidelberg, Berlin, Heidelberg. Calvete, H. I., Gal´e, C. and Iranzo, J. A. (2016) “Meals: A multiobjective evolutionary algorithm with local search for solving the bi-objective ring star problem”. European Journal of Operational Research, vol. 250(2): 377–388. Ghiani, G. and Improta, G. (2000) “An efficient transformation of the generalized vehicle routing problem”. European Journal of Operational Research, vol. 122(1): 11–17. Haimes, Y. (1971) “On a bicriterion formulation of the problems of integrated system identification and system optimization”. IEEE transactions on systems, man, and cybernetics, (3): 296–297. Hoshino, E. A. and De Souza, C. C. (2012) “A branch-and-cut-and-price approach for the capacitated m-ring–star problem”. Discrete Applied Mathematics, vol. 160(18): 2728–2741. Jolfaei, A. A., Alinaghian, M., Bahrami, R. and Tirkolaee, E. B. (2023) “Generalized vehicle routing problem: Contemporary trends and research directions”. Heliyon. Labb´e, M., Laporte, G., Mart´ın, I. R. and Gonz´alez, J. J. S. (2004) “The ring star problem: Polyhedral analysis and exact algorithm”. Networks: An International Journal, vol. 43(3): 177–189. Lana, I. V. (2016) “El problema del m-anillo estrella”. Tesis de la Universidad de Zaragoza. Leite, B. S. C. (2023) “The vehicle routing problem: Exact and heuristic solutions”. https://towardsdatascience.com/ the-vehicle-routing-problem-exact-and-heuristic-solutions-c411c0f4d734. Naji-Azimi, Z., Salari, M. and Toth, P. (2012) “An integer linear programming based heuristic for the capacitated m-ring-star problem”. European Journal of Operational Research, vol. 217(1): 17–25. Oakey, A., Martinez-Sykora, A. and Cherrett, T. (2023) “Improving the efficiency of patient diagnostic specimen collection with the aid of a multi-modal routing algorithm”. Computers & Operations Research, vol. 157: 106265. 55 P´erez, J. A. M., Moreno-Vega, J. M. and Martın, I. R. (2003) “Variable neighborhood tabu search and its application to the median cycle problem”. European Journal of Operational Research, vol. 151(2): 365–378. Pozo, M. A., Puerto, J. and Rold´an, I. (2024) “The biobjective minimum-cost perfect matching problem and chinese postman problem”. International Transactions in Operational Research, vol. 31(5): 2988–3009. Przybylski, A., Gandibleux, X. and Ehrgott, M. (2010) “A recursive algorithm for finding all nondominated extreme points in the outcome set of a multiobjective integer programme”. INFORMS Journal on Computing, vol. 22(3): 371–386. P¨usk¨ul, O. M., Aslan, D., Onay, C., Erdogan, M. S. and Ta¸sgetiren, M. F. (2022) “Designing a railway network in cesme, izmir with bi-objective ring star problem”. In N. M. Durakbasa and M. G. Gen¸cyılmaz, eds., Digitizing Production Systems, pp. 665–674. Springer International Publishing, Cham. Simonetti, L., Frota, Y. and de Souza, C. C. (2011) “The ring-star problem: a new integer programming formulation and a branch-and-cut algorithm”. Discrete applied mathematics, vol. 159(16): 1901– 1914. Toth, P. and Vigo, D. (2002) The Vehicle Routing Problem. SIAM. ISBN 0-89871-498-2. Wolsey, L. A. (2020) Integer programming. John Wiley & Sons. ISBN 9781119606536. 56