Câmpus de São José do Rio Preto Eduardo dos Santos Teixeira Formulations and a heuristic approach for logistics problems with d-relaxed priority rule São José do Rio Preto 2024 Eduardo dos Santos Teixeira Formulations and a heuristic approach for logistics problems with d-relaxed priority rule Tese apresentada como parte dos requisitos para obtenção do título de Doutor em Matemática, junto ao Programa de Pós-Graduação em Matemática, do Instituto de Biociências, Letras e Ciências Ex- atas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto. Financiadora: CAPES Orientador: Prof. Dr. Silvio Alexandre de Araujo São José do Rio Preto 2024 T266f Teixeira, Eduardo dos Santos Formulations and a heuristic approach for logistics problems with d-relaxed priority rule / Eduardo dos Santos Teixeira. -- São José do Rio Preto, 2024 117 p. : il., tabs. Tese (doutorado) - Universidade Estadual Paulista (Unesp), Instituto de Biociências Letras e Ciências Exatas, São José do Rio Preto Orientador: Silvio Alexandre de Araujo 1. Programação Inteira Mista. 2. Problema do Caixeiro Viajante. 3. Problema de Roteamento de Veículos. 4. Prioridades. 5. Heurística. I. Título. Sistema de geração automática de fichas catalográficas da Unesp. Biblioteca do Instituto de Biociências Letras e Ciências Exatas, São José do Rio Preto. Dados fornecidos pelo autor(a). Essa ficha não pode ser modificada. Eduardo dos Santos Teixeira Formulations and a heuristic approach for logistics problems with d-relaxed priority rule Tese apresentada como parte dos requi- sitos para obtenção do título de Doutor em Matemática, junto ao Programa de Pós- Graduação em Matemática, do Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto. Financiadora: CAPES Comissão Examinadora Prof. Dr. Silvio Alexandre de Araujo Orientador Profa. Dra. Maria do Socorro Nogueira Rangel UNESP - Câmpus de São José do Rio Preto Prof. Dr. Pedro Augusto Munari Junior UFSCAR - Câmpus de São Carlos Profa. Dra. Maria José Pinto IEAv/CTA - São José dos Campos Prof. Dr. Alf Kimms Universität Duisburg-Essen - Duisburgo, Alemanha São José do Rio Preto 06 de março de 2024 À minha família, dedico AGRADECIMENTOS A presente tese de doutorado é fruto de uma trajetória de 12 anos de Universi- dade, os quais são, por sua vez, a sequência de outros tantos anos de estudos e vivências que se deram antes. Dessa forma, correndo o risco de, injustamente, me esquecer de algumas pessoas que foram parte importante desta conquista, gostaria de registrar aqui os meus sinceros agradecimentos a todos aqueles que estiveram próximos, de qualquer forma, durante o desenvolvimento deste trabalho. Agradeço, em primeiro lugar, a Deus, inteligência maior e causa primária de to- das as coisas, pela oportunidade de chegar até aqui, pela coragem e pelo ânimo de seguir adiante com o trabalho, mesmos nos dias de cansaço, aflição ou ansiedade que insistiram em acontecer. Agradeço, também, a toda a minha família, que sempre me apoiou em minha trajetória acadêmica e sempre me ensinou que a educação é uma chave com o poder de abrir grandes portas na vida. Através de palavras e exemplos valiosos, essa lição me motivou a buscar subir os degraus da vida acadêmica e a trabalhar pela conquista deste doutorado. Em particular, agradeço aos meus pais, Rita e Matias, ao meu irmão Augusto, e aos meus avós Alice a Alcides. Agradeço à Ieda, minha namorada, por todo o apoio, toda a comprensão e todo o incentivo durante a realização deste trabalho, bem como de todos os outros planos e projetos que se realizaram durante o desenvolvimento desta tese. Agradeço ao Programa de Pós Graduação em Matemática (PPGMAT) e aos professores dos Departamentos de Matemática (DMAT) e de Matemática Aplicada (DMAP) do Ibilce - em especial ao meu orientador, professor Sílvio, cujo trabalho foi fundamental para que cada peça deste estudo fosse confeccionada e colocada no seu devido lugar. Sem a sua coordenação e os seus apontamentos, com certeza esta tese perderia em organização, clareza e conteúdo. Agradeço, também, aos demais funcionários do Ibilce, sem os quais nenhuma atividade da universidade poderia se desenvolver de forma satisfatória. Em parti- cular, deixo o meu agradecimento ao Thiago, que me ajudou durante esses anos com o acesso aos recursos de infra-estrutura utilizados nos testes computacionais realizados. Quero agradecer também a todos os membros da banca de defesa deste tra- balho, os quais aceitaram o convite para avaliar essa tese e colaborar com a sua elaboração final. Dentre uma lista vasta de outras pessoas - amigos, colegas e parceiros - que estiveram ao meu lado durante a realização deste curso, parcial ou totalmente, ou que me ajudaram de alguma forma a chegar até aqui, deixo aqui o meu agrade- cimento à Letícia, à Gi, à Ju, à Ana, à Mariele, ao Luis, ao Felipe, ao Décio e a tantos outros que me apoiaram nessa trajetória. O tópico central desta tese foi proposto e inicialmente desenvolvido como fruto de um trabalho em parceria com as seguintes pessoas - Drielly, Natália e Marcos - às quais deixo um agradecimento em particular. O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoa- mento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001. O homem que julga infalível a sua razão está bem perto do erro, Allan Kardec RESUMO Esta tese considera extensões de dois problemas clássicos de pesquisa opera- cional no campo da logística em que regras de prioridade são incorporadas às versões clássicas, dando origem ao Problema do Caixeiro Viajante Clusterizado com Regra de Prioridade d-Relaxada (CTSP-d) e ao Problema de Roteamento de Veículos Clusterizado com Regra de Prioridade d-Relaxada (CluVRP-d), que corres- ponde a um caso particular do Problema de Roteamento de Veículos com Regras de Prioridade Relaxadas (VRP-RPR). Em ambos os problemas, a ordem em que os vértices são visitados é relevante porque é necessário considerar, além da dis- tância percorrida, algum tipo de prioridade de visita entre os vértices, com base em diferentes situações reais. Neste estudo, formulações propostas na literatura para o CTSP-d são melhoradas utilizando-se desigualdades válidas, bem como novas for- mulações baseadas em variáveis de precedência são propostas. Outra contribuição desta tese é expandir a literatura sobre o VRP-RPR propondo o CluVRP-d, uma ver- são do problema onde todos os vértices devem ser visitados exatamente uma vez e os veículos possuem limitações de capacidade. Duas novas formulações matemáti- cas são propostas para o problema, considerando os casos Local e Global Timing, e uma abordagem heurística é apresentada para tratar de instâncias maiores do caso Local Timing. Por fim, resultados computacionais, baseados em dados da literatura, são apresentados para demonstrar a competitividade das formulações propostas quando resolvidas por um pacote de otimização e, no caso do CluVRP-d, comparamos a formulação do caso Local Timing com os resultados obtidos através da abordagem heurística. Palavras-chave: Programação Inteira Mista, Problema do Caixeiro Viajante, Prob- lema de Roteamento de Veículos, Prioridades, Heurística. ABSTRACT This thesis considers extensions of two classical operations research problems in the logistics field where priority rules are embedded to the classical versions, defin- ing the Clustered Traveling Salesman Problem with d-Relaxed Priority Rule (CTSP- d) and the Clustered Vehicle Routing Problem with d-Relaxed Priority Rule (CluVRP- d), that corresponds to a particular case of the Vehicle Routing Problem with Re- laxed Priority Rules (VRP-RPR). In both problems, the order in which the vertices are visited is relevant because it is necessary to consider, in addition to the dis- tance traveled, some kind of visiting priorities among the vertices, based on different real situations. In this study, formulations for the CTSP-d are improved using valid inequalities as well as new formulations based on precedence variables are pro- posed. Another contribution of this thesis is to extend the literature on VRP-RPR by proposing the CluVRP-d, a version of the problem where all vertices must be visited exactly once and the vehicles have capacity limitations. Two new mathematical for- mulations are proposed considering local and global timing cases, and a heuristic approach is presented to tackle bigger instances for the local timing case. Finally, computational results, based on data from the literature, are presented to demon- strate the competitiveness of the proposed formulations running the models on an optimization package and, in the CluVRP-d case, we compare the Local Timing case formulation to the results obtained by the heuristic approach. Keywords: Mixed Integer Programming, Traveling Salesman Problem, Vehicle Rout- ing Problem, Priorities, Heuristic. List of Figures 5.1 H2020, MTZ1, GP1, SSB1 and SST1 runtime performance profile for smaller instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 H2020, MTZ1, GP1, SSB1 and SST1 GAP performance profile for smaller instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.3 H2020, MTZ2, GP2, SSB2 and SST2 runtime performance profile for smaller instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.4 H2020, MTZ2, GP2, SSB2 and SST2 GAP performance profile for smaller instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.5 H2020, MTZ2, GP2, SSB2 and SST2 runtime performance profile for 100 vertices instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.6 H2020, MTZ2, GP2, SSB2 and SST2 GAP performance profile for 100 vertices instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.7 The berlin52 TSP instance and its solution . . . . . . . . . . . . . . . . 76 5.8 Optimal tours obtained for different values of d for the berlin52-R-5-d-b instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.9 Line plots showing how the optimal tour length decreases as the value of d increases for some selected instances . . . . . . . . . . . . . . . . . 78 List of Tables 3.1 Size order and reference numbers of each formulation . . . . . . . . . . 42 5.1 Average values of Objective, Runtime and GAP, and number of Best Known Solutions (BKS) found by the solver for the smaller instances running formulations H2020, MTZ1, GP1, SSB1 and SST1 . . . . . . . 65 5.2 Average values of Objective, Runtime and GAP, and number of Best Known Solutions (BKS) found by the solver for the smaller instances running formulations H2020, MTZ2, GP2, SSB2 and SST2 . . . . . . . 68 5.3 Average values of Objective, Runtime and GAP, and number of Best Known Solutions (BKS) found by the solver for the 100 areas instances running formulations H2020, MTZ2, GP2, SSB2 and SST2 . . . . . . . 71 5.4 Solver results from the H2020 and MTZ2 formulations for the 200-city random instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.5 Solver results from the H2020 and MTZ2 formulations for the 200-city clustered instances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.6 Fulfilling priorities in the selected instances . . . . . . . . . . . . . . . . 79 6.1 Summary of results for the MILP models proposed, showing the means of objective value, runtime and GAP for each vehicle fleet. . . . . . . . 85 6.2 Average objective value, average runtime and percentage of instances the best solution was found for each solution approach on the smaller instances considering the Local Timing Model. . . . . . . . . . . . . . . 86 6.3 Average best objective value and total runtime when using the heuristic approach on the instances with 100 and 200 vertices. . . . . . . . . . . 87 A.1 Solver results from the H2020, MTZ1, GP1, SSB1 and SST1 formula- tions for the smaller random instances . . . . . . . . . . . . . . . . . . 101 A.2 Solver results from the H2020, MTZ1, GP1, SSB1 and SST1 formula- tions for the smaller clustered instances . . . . . . . . . . . . . . . . . . 102 A.3 Solver results from the H2020, MTZ2, GP2, SSB2 and SST2 formula- tions for the smaller random instances . . . . . . . . . . . . . . . . . . 103 A.4 Solver results from the H2020, MTZ2, GP2, SSB2 and SST2 formula- tions for the smaller clustered instances . . . . . . . . . . . . . . . . . . 104 A.5 Solver results from the H2020, MTZ2, GP2, SSB2 and SST2 formula- tions for the 100-city random instances . . . . . . . . . . . . . . . . . . 105 A.6 Solver results from the H2020, MTZ2, GP2, SSB2 and SST2 formula- tions for the 100-city clustered instances . . . . . . . . . . . . . . . . . 106 A.7 Small random instances - Local Timing Model Case . . . . . . . . . . . 108 A.8 Small clustered instances - Local Timing Model Case . . . . . . . . . . 109 A.9 Small random instances - Global Timing Model Case . . . . . . . . . . 110 A.10 Small clustered instances - Global Timing Model Case . . . . . . . . . 111 A.11 Small random instances - Local Timing Model Heuristic . . . . . . . . . 112 A.12 Small cluster instances - Local Timing Model Heuristic . . . . . . . . . 113 A.13 100 vertices instances - Local Timing Model Heuristic . . . . . . . . . . 114 A.14 200 vertices instances - Local Timing Model Heuristic . . . . . . . . . . 115 List of Abbreviations CTSP-d Clustered Traveling Salesman Problem with d-Relaxed Priority Rule CluVRP-d Clustered Vehicle Routing Problem with d-Relaxed Priority Rule MILP Mixed-integer linear programming TSP Traveling Salesman Problem VRP-RPR Vehicle Routing Problem with Relaxed Priority Rules VRP Vehicle Routing Problem Contents 1 Introduction 15 2 Literature review 18 2.1 Clustered Traveling Salesman Problem . . . . . . . . . . . . . . . . . . 18 2.2 Vehicle Routing Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.1 VRPs with priorities . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.2 VRP with clustering . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3 Formulations for the CTSP-d 25 3.1 Miller-Tucker-Zemlin based formulations . . . . . . . . . . . . . . . . . 27 3.1.1 Proposed valid inequalities for the MTZ1 formulation . . . . . . 28 3.2 Precedence variable-based formulations . . . . . . . . . . . . . . . . . . 33 3.2.1 Proposed valid inequalities for the precedence variables-based formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4 Formulations and a heuristic approach for the CluVRP-d 43 4.1 Formulations for the CluVRP-d . . . . . . . . . . . . . . . . . . . . . . 45 4.1.1 Local Timing Model . . . . . . . . . . . . . . . . . . . . . . . . 46 4.1.2 Global Timing Model . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1.3 Valid Inequalities . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2 A heuristic approach for the Local Timing case . . . . . . . . . . . . . 49 4.2.1 Clustering Phase . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2.2 Initial Feasible Solution Phase . . . . . . . . . . . . . . . . . . . 53 4.2.3 Improvement Phase . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2.4 The complete heuristic procedure . . . . . . . . . . . . . . . . . 59 5 Computational Results for the CTSP-d 61 5.1 Benchmark instances and experimental settings . . . . . . . . . . . . . 61 5.2 Comparison of the results from the formulations H2020, MTZ1, GP1, SSB1 and SST1 for smaller instances . . . . . . . . . . . . . . . . . . . 63 5.3 Comparison of the results from the formulations H2020, MTZ2, GP2, SSB2 and SST2 for smaller instances . . . . . . . . . . . . . . . . . . . 67 5.4 Comparison of the results from the formulations H2020, MTZ2, GP2, SSB2 and SST2 for large instances . . . . . . . . . . . . . . . . . . . . 70 5.5 How the TSP, CTSP-PO and the CTSP-d are related . . . . . . . . . . 75 6 Computational Results for the CluVRP-d 82 6.1 Benchmark instances and experimental settings . . . . . . . . . . . . . 82 6.2 Computational results obtained by the proposed formulations . . . . . 84 6.3 Computational results obtained by the proposed heuristic . . . . . . . 86 7 Conclusion and future research 88 References 91 Appendix A Numerical results 100 A.1 CTSP-d Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 A.2 CluVRP-d Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 1 Introduction In recent decades, logistics related problems have been studied in numerous papers, with different mathematical modeling and solution approaches. A central problem in this context is the Traveling Salesman Problem (TSP), also called the Traveling Sales- person Problem - one of the most important problems in combinatorial optimization and computer science (Dantzig et al., 1954). In its classical version, the TSP may be described in the following way: given a set of cities, and knowing the distance between each of these cities, we need to find the shortest possible route such that each city is visited only once, including the city of origin where the route also ends. In other words, the objective is to find the round trip with the minimum possible length. This is known to be an NP-hard problem (Garey and Johnson, 1979) and so, to solve it to optimality, finding and proving that a solution is the best one is not always possible in practice. In such cases, the objective of the problem may be to build the shortest possible tour within computational or solution-time restrictions. An important gen- eralization of this problem is the Vehicle Routing Problem (VRP), where instead of considering only one traveling salesperson visiting a set of cities, the VRP is defined by a fleet of vehicles, and so the problem is to find the optimal routes for the fleet of vehicles, delivering goods or services to a set of customers while minimizing the total cost or distance traveled. A common and useful way to model the TSP and the VRP is by using an oriented graph, where all the depot from where the vehicles leave and all the customers that must be visited are represented by vertices, and the possibility of traveling among two areas is represented by an arc. In the classical version of the problems we consider only one starting depot and that is possible to travel from any area (depot or customer) to any different area, what leads to a complete graph. In this approach, the problems can be seen as the challenge of choosing the best set of arcs that satisfy all constraints. In this thesis we will consider this modeling resource for the mathematical formulations proposed for all problems, and so we may use the words vertex and arcs to describe the problems from here on. On the classical TSP and VRP definitions, the order in which the customers are visited has no restrictions. In some real situations, just this condition may not be 15 16 Introduction sufficient to represent the problem, as there are cases in which some kind of visiting priorities among the customers must be considered - in general, customers with the highest priorities need to be served before lower priority ones. According to Doan et al. (2021) this problem has several applications where different priority groups are needed, such as humanitarian relief operations (Panchamgam, 2011), unmanned aerial vehicles (Shetty et al., 2008), distribution of commercial products when the priority is determined by the customers’ storage level and technician scheduling problems faced by internet service providers (Hà et al., 2020). A deeper discussion on variants of the TSP and VRP involving priorities is developed in Chapter 2. The focus of the present thesis is to study a prioritization rule called “d-relaxed priority rule” in the contexts of the TSP and the VRP separately. To apply the rule, first we need to assign priorities to each vertex of the problem in descending order, that is, vertices that receive priority 0 have the highest priority level, followed by vertices with priority 1, priority 2, and so on. The d-relaxed priority rule is used to manage the operational flexibility when visiting the vertex by allowing the vehicles to visit a whole set of similar priority level sets in an interchangeable manner but forbidding visits to lower priority sets until the higher priority ones have been totally served. By setting d = 0, all the priorities must be strictly verified, which usually leads to costly solutions, while higher values of d allows the solution to ignore partially (or even totally) the priorities, leading to smaller/cheaper routes. We will better describe the details of the rule in Chapter 3 for the TSP and Chapter 4 for the VRP. Two important and still not deeply explored logistics problems with priorities were described in Panchamgam (2011), called Hierarchical Traveling Salesman Problem with d-relaxed priority rule (HTSP-d) and the Vehicle Routing Problem with Relaxed Pri- ority Rules (VRP-RPR). Doan et al. (2021) study VRP-RPR in the context of two important prioritization rules, namely, the d-relaxed priority and the Order of Demand Fulfilment (ODF), where the second rule is used to ensure that at the end of the plan- ning horizon, all the higher priority vertices need to be visited before the vehicles can visit any lower priority one. The d-relaxed priority rule was also studied as a prioriti- zation rule for the TSP (Panchamgam, 2011; Hà et al., 2020; Teixeira and de Araujo, 2024). In this context, since all the vertices need to be visited, we do not need to consider the ODF rule, and so the problem reduces to a TSP with d-relaxed priority rule. This problem is called HTSP-d or the Clustered Traveling Salesman problem with d-relaxed priority rule (CTSP-d). In this thesis we propose a generalization of the CTSP-d as a problem with a set of vehicles, defining the Clustered Vehicle Routing Problem with d-relaxed Priority Rule (CluVRP-d), which is a particular case of the VRP-RPR. Introduction 17 The research in this thesis is inspired by the agricultural pest control presented in Teixeira et al., 2021, where the authors describe an ant control problem in forestry that must be handle considering priorities issues. The problem can be described in the following way: given a set of forest areas, a working group must visit regularly each one of these areas to develop a pest control service, combating an active infestation. Depending on the severity of the infestation, it may impair the forest production in different levels, and so the working group must visit all areas considering the individual urgency levels, which are grouped into priority visiting sets based on technical analysts criteria. This problem is identified with the CTSP-d problem, and the multi working group version of the problem can be seen as the CluVRP-d. The remainder of the text is organized in the following way: in Chapter 2 we present a literature review on the topic logistics problems with priority constraints, discussing different problems proposed involving the TSP and the VRP in this context. Chapter 3 presents formulations for the CTSP-d based in the classical Miller et al. (1960) formulation and in precedence variables, showing different MILP models for the problem and related results. In Chapter 4 are introduced formulations for the CluVRP- d considering the two ways the d-relaxed priority rule can be applied in the VRP - namely, the Local Timing Model and the Global Timing Model. Valid inequalities are presented for the models, and a heuristic approach for the Local Timing case is proposed. Chapters 5 and 6 present, respectively, the computational results obtained running the formulations proposed for the CTSP-d in a commercial MILP solver, and the computational results obtained for the CluVRP-d running both the formulations using a MILP solver and the heuristic proposed for the problem. Finally, in Chapter 7 are the conclusions and future research proposals. It is worth mentioning that a version of the contents of Chapters 3 and 5 was published on the paper Teixeira and de Araujo (2024). 2 Literature review This literature review presents a brief overview of the problems and papers related to our research. The review is divided into applications of prioritization rules or related topics to the TSP and the VRP, separately. This is not supposed to be an exhaustive review about TSP/VRP variants with priority criteria, but a list of the recent works related at some level to them. 2.1 Clustered Traveling Salesman Problem An important problem that has been widely studied in the literature is the so-called Clustered Traveling Salesman Problem (CTSP) introduced by Chisman (1975). In this problem, the traveling salesman not only must visit each vertex exactly once, but the vertices are partitioned into disjoint subsets, called clusters, where all vertices included in a cluster must be visited contiguously. In the CTSP, there is no order relationship between the vertices in each cluster, as well as no order in which the clusters must be visited. In Laporte et al. (1997) and Potvin and Guertin (1998) the authors consider the CTSP with a prespecified order on the clusters and introduced the Clustered Traveling Salesman Problem with a Prespecified Order (CTSP-PO). In this problem, there are the sets of vertices with associated priorities defined by V0, V1, ..., VPmax , but no d parameter. Since this problem has a prespecified order on the clusters visiting that must be strictly followed, any feasible solution for the CTSP-PO must start at the origin vertex 0 followed by some vertex in V0, which is the highest priority set. Vertices in V0 must be visited contiguously, and once the last one is left, the next vertex must be in V1, which is the next highest priority set still not visited. This process repeats until all vertices in VPmax have been visited. We see a hard priority fulfillment on this behavior, that may be relaxed in some practical situations. A natural application of the CTSP-PO occurs when vertices requiring the same level of urgency can be clustered together, and higher priority clusters should be visited before others. Several authors have developed solution approaches to the CTSP (see, for example, Mestria (2018)). The CTSP-d is defined in a similar but more general way than the CTSP-PO, since 18 Clustered Traveling Salesman Problem 19 there is flexibility in the prespecified order the clusters must be served, which is defined by the d parameter of the rule. The CTSP-PO correspond to the CTSP-d when we take d = 0 and there is no flexibility allowed on the clusters visiting order. Another related problem is the Traveling Salesman Problem with Sequence Priori- ties (TSPwSP), introduced by Schmitz and Niemann (2009). The TSPwSP is defined by assigning priorities to each vertex and measuring the violations of these priorities using a penalty function. The defined priorities stand for the order in which each vertex asked for the salesman visit, which in practice may represent a certain service that must be provided to a set of clients, adding a “first-come-first-served” format to the problem. Since the objectives of minimizing the total cost of the tour and mini- mizing the violations on the service order line are conflicting in general, the authors proposed a bicriterion approach to manage the trade-off between them. When compar- ing the TSPwSP and the CTSP-d, there are clearly some important differences. First, the TSPwSP defines a bi-objective model, in which the priorities are assigned to the vertex, to define the “first-come-first-served” policy intended by the authors while, in the CTSP-d, the only objective is minimizing the total distance of the tour while the priorities are considered as constraints over the whole clusters of vertices. Moreover, the TSPwSP can ignore part of the priorities, since the model only deals with allowing the vertices to be more or less far from the “ideal positions” they should be. Similarly, in the Travelling Salesman Problem with Priority Prizes (TSPPP), de- scribed in Pureza et al. (2018), instead of adding a penalty objective function to the TSP, the authors propose a model based on the Quadratic Assignment Problem, in which they include a prize according to the priorities in the model, which may be collected by the salesman by meeting these order visiting prizes. In the Multicommod- ity Traveling Salesman Problem (MTSP), described by da Silva et al. (2019), when connecting two vertices (customers), product types and quantities being transported are considered in the objective function and priorities are defined based on the size of demands or on the value/risk of the products that will be delivered. Once again, when compared to the CTSP-d, both problems (TSPPP and MTSP) are not imposing pri- orities that must be respected in some sense, but just adding conditions on preferable positions some vertices should assume in the tour, and so the solutions obtained via TSPPP or MTSP may not be the ideal for situations where some vertices need to be visited before others. More specifically, when dealing with humanitarian logistics, in planning relief operations after a natural disaster, the objective is clearly not to obtain financial benefits but to save people’s lives. There are other important TSP variants, such as the Generalized Traveling Sales- man Problem (GTSP), discussed in Srivastava et al. (1969) and modeled in Laporte 20 Literature review and Nobert (1983), Laporte et al. (1987) and Pop (2007), in which the vertices are partitioned into clusters, as in the CTSP but, in the GTSP, the objective is to build a tour in which the salesman visits only one vertex in each cluster. There are exact and heuristic methods to solve this problem (see, for example, Pintea et al. (2017) and Cosma et al. (2021)), and there is a clustered version of this problem, namely the Clustered Generalized Traveling Salesman Problem (CGTSP - see Baniasadi et al. (2020)). An extension of the GTSP is the Family Traveling Salesman Problem (FTSP), presented in Bernardino and Paias (2018), where the vertices are divided into clusters, called families in the context of the problem, and the objective is to determine a mini- mum cost route in which the traveling salesman visits a given number of cities in each family. Another relevant way the vertices in the TSP can be ordered is the Traveling Salesman Problem with Time Windows (TSPTW) (Dumas et al., 1995), that does not consider priorities among the vertices explicitly but the time windows in which each customer must be visited can represent the priority among the vertices. In practice, the TSPTW imposes ranges of positions each vertex may figure in a solution, but this is not sufficient to define the d-relaxed priority rule. It is also worth mentioning the TSP variant so-called Crew Scheduling and Routing Problem (CSRP) (Maya Duque et al., 2016; Moreno et al., 2019, 2020) that consists of finding a single route (TSP) to visit nodes that need to be repaired after a disaster and the order of visits is important and implied by the need of relief paths reaching the damaged areas as soon as possible. The Precedence Constrained Traveling Salesman Problem (PCTSP), also called the Sequential Ordering Problem (Escudero, 1988), is another relevant TSP variant, where there are precedence relationships between part of the vertices of the graph of the problem. It introduces a new set of constraints in the problem, imposing that the traveling salesman visit some defined vertices before (but not necessarily immediately before) other ones. The PCTSP can be extended into the Precedence Constrained Generalized Traveling Salesman Problem (PCGTSP), which brings together both the conditions of GTSP and PCTSP, defining a problem where the vertices are divided into clusters, and the traveling salesman must visit each cluster in a single vertex, respecting the precedence conditions imposed (Salman et al., 2020; Khachai et al., 2023). Some recent optimal solutions and best bounds for the PCTSP can be find in Cire and van Hoeve (2013) and Gouveia and Ruthmair (2015), and some recent advances in formulations for the problem are presented in Letchford and Salazar González (2016) and Gouveia et al. (2018). In Chapter 3 we discuss how the CTSP-d can be seen as a particular case of the PCTSP and how to adapt some classical formulations for it. This can be an important step to allow the use of other models and efficient solution methods for the PCTSP as modeling/solving resources for the CTSP-d. Vehicle Routing Problem 21 In this study, new formulations for the CTSP-d are presented that, as discussed above, are different from TSP with prize collection models or models imposing fixed or desirable positions for each vertex. They can be seen as a natural generalization of the CTSP-PO for problems which consider trade-offs between tour costs and the hardness of the priority assigned to the clusters. Although the consideration of priorities is a relevant topic in the routing literature, to the best of our knowledge, the d-relaxed priority rule has not been totally explored. As mentioned before, this rule was proposed by Panchamgam (2011) in the context of retail operations and humanitarian logistics, and Hà et al. (2020) proposed a reformulation of the original model, leading to a more precise representation of the priority rule and better computational results. Teixeira and de Araujo (2024) presented new formulations for the problem based in the work of Hà et al. (2020), using the classical Miller et al. (1960) approach to model the problem, as well as precedence variables models. Doan et al. (2023) extended the work of Hà et al. (2020) by proposing alternative formulations and an iterated local search. 2.2 Vehicle Routing Problem Considering the VRP, the vast literature has many different ways of treating the existence of priorities on the solutions, or any other kind of vertex ordering requirements that use some criteria to distinguish transportation service between customers such as: priority assignment, weighting factor or penalty objective. On the other hand, the VRP with d-relaxed priority rule was first described in Panchamgam (2011) in the context of disaster relief operations as a way to balance the trade-off between the priorities fulfillment and the costs involved in building the routes in the VRP. To the best of our knowledge, the only paper dealing with the VRP-RPR is Doan et al. (2021), where the authors solved it heuristically. The goal of this literature review is, in this sense, to show some of the strategies used to treat priorities in the VRP literature, comparing it with the d-relaxed priority rule, and to analyze how these approaches differ from each other. To do so, we first present other variants of the VRP with priorities in Subsection 2.2.1; next, we consider VPR with clustering approaches in Subsection 2.2.2 and, finally, some application that consider priorities in Subsection 2.2.3. 2.2.1 VRPs with priorities Several extensions of the VRP are related to priorities and have been studied in the literature, such as: the VRP with service levels (VRP-SL) (Bulhões et al., 2018), the Inventory Routing Problem (IRP) (Archetti and Ljubić, 2022) and the VRP with Profits (VRPPs) (Archetti et al., 2014). The VRP-SL takes the viewpoint of a third 22 Literature review party logistics provider and because of limited fleet size and vehicle capacity, it is necessary to select a subset of deliveries and determine cost-efficient vehicle routes in such a way that the overall activity is profitable and that the agreements with the partners are respected. In the VRPPs, the set of customers to serve is not given and needs to be decided according to its profit, as well as, how to cluster the customers to be served in different routes. Thus, the routes can be measured both in terms of cost and in terms of profit. According to Doan et al. (2021), when demand is replaced by profit, the VRP-RPR can be reduced to some important problems, such as the Team Orienteering Problem (TOP) (Chao et al., 1996; Boussier et al., 2007), that is a special case of the VRP, where there are a set of potential customers to be visited and a fleet of capacitated vehicles to be assigned to them. Due to some practical limitations, it is considered that not all the customers can be visited during the tours built, and so the objective is to construct vehicle routes that maximize the reward obtained by visiting each vertex. Other relevant related extensions are: the Dynamic VRP With Time Windows and Multiple Priorities (Yang et al., 2015) and the Multi-objective Vehicle Routing Problem with Multiple Prioritized Time Windows (Beheshti et al., 2015). In the former, there is a fleet of vehicles that is used to serve customers with time windows and different priority levels. Customers with higher priority should have a higher quality of service. Additionally, part of customers are revealed dynamically during the execution of the routes. The latter article, in addition to the minimization of the total traveling cost, aims to maximize customer satisfaction regarding prioritized time windows. 2.2.2 VRP with clustering Some VRP variants consider that the customers are also divided into clusters and are served according to the cluster to which they belong. These problems are related to those with priorities because in several applications the cluster can be defined according to the priority of each customer and, in fact, we will use this strategy in Chapters 3 and 4. The Generalized Vehicle Routing Problem (GVRP) (Ghiani and Improta, 2000; Hà et al., 2014), similarly to the GTSP, considers that each cluster contains a number of vertices, and a tour includes exactly one vertex from each cluster and satisfies capacity constraints. The Selective Vehicle Routing Problem (SVRP) (Sabo et al., 2020) is a generalization of the GVRP where customers may belong to one or more clusters. The Clustered Vehicle Routing Problem (CluVRP), studied in Battarra et al. (2014) and Pop et al. (2018), is the multi-vehicle case of the CTSP, where once a vehicle visits one customer in a cluster it must visit all the remaining customers of that cluster before leaving it. In Expósito-Izquierdo et al. (2016) the authors consider the capacitated Vehicle Routing Problem 23 version on the CluVRP. Yahiaoui et al. (2019) present the Clustered Team Orienteering Problem (CluTOP) where a profit is associated with each cluster and is gained only if all of its customers are served considering a set of available vehicles with a limited travel time. It is also important to highlight some papers that consider clustering in a more general way, meaning assigning similar objects to groups. These problems have been considered in a broad range of applications and, in several cases, the size of the clusters is limited; consequently these capacitated clustering problems have received increasing attention (Stefanello et al., 2015; Martínez-Gavara et al., 2017; Chaves et al., 2018; Zhou et al., 2019; Gnägi and Baumann, 2021). 2.2.3 Applications One of the main sources of application of VRP with priorities rules are humanitar- ian logistics (Campbell et al., 2008) including emergency response (Bretschneider and Kimms, 2011; Jacobson et al., 2012). For a complete analysis about the differences between commercial and humanitarian logistics the reader can refer to Beamon and Balcik (2008), Gralla and Goentzel (2018) and Wassenhove (2006). The literature on logistics for disaster relief operations deals with the existence of priorities in different ways, and this is a very relevant topic, since in this kind of application, bad planning may intensify the suffering of the affected people, exposing them to injuries, diseases or even the risk of death. Barbarosoglu et al. (2002) present two models to deal with tactical and operational planning in the context of routing a fleet of helicopters to work on disaster relief operations. Both the proposed models are mixed integer and the authors describe a hierarchical multi-criteria approach to relate the two stages of the planning. The objective is to minimize the total tour duration for the helicopters, and no explicit priority constraints are defined in any of the models. A problem similar to the VRP is described by Bish (2011), who models the so-called “Bus Evacuation Problem” (BEP), in which the vertices are divided in three subsets: yard vertices, where the buses are initially located and do not return during their tours, demand vertices, where the buses need to visit to do the evacuation work, and shelter vertices, the places where people are taken after being rescued. The problem is similar to a VRP with multi-depot and split delivery, and the authors consider two distinct objectives: to minimize the maximum tour cost over the buses, and to minimize the total cost of the routing problem. The solution procedure proposed is a two-phase heuristic to build and improve a solution using ideas similar to the ones commonly used in the context of the VRP. Briskorn et al. (2020) consider an integrated problem that consists in unblocking roads in order to make demand locations accessible and delivering relief 24 Literature review goods in order to satisfy demand considering strict deadlines. Other applications in the disaster relief context including emergency response can be found in Altay and Green (2006), Sheu (2007), Sheu (2010), Chiu and Zheng (2007), Özdamar and Demir (2012), Simpson and Hancock (2009) and Kimms and Maiwald (2018). 3 Formulations for the CTSP-d In this chapter we will present formulations for the CTSP-d based in the classical Miller et al. (1960) formulation and in precedence variables. Before presenting the formulations we will formally define the d-relaxed priority rule considering the following sets and parameter. • Sets: V = {0, 1, 2, ..., n − 1}: Set of the n vertices of the problem, where the vertex 0 represents the start city of the tour, and the other vertices represent the cities that must be visited; A ⊆ {(i, j) : i, j ∈ V, i ̸= j}: Set of the arcs of the problem, representing the possibility of the traveling salesman to travel from vertex i to vertex j. For the sake of simplicity, in this thesis, we assume a complete graph, meaning that the set A contains arcs linking all the different vertices in the problem. P = {0, 1, 2, ..., Pmax}: Set of the different priority levels of the vertices in the problem; Vp: Set of vertices from V \{0} with priority level p ∈ P . The sets Vp define a partition in V \{0}, i.e., the sets are disjoints and the union of them all is equal to V \{0}, and are ordered in descending order with respect to the priority level they represent, that is, set V0 include the vertices with the highest priority level, followed by V1 up to VPmax . • Parameter: d: Parameter used to define the hierarchy relaxation level over the priority sets. For the d-relaxed priority rule, the first vertex after the origin may be in any set among V0, ..., Vd, and the vertices on these sets may be visited in an interchangeable way, disregarding the different priority levels. Vertices in sets Vd+1, ..., VPmax stay forbidden until all vertices in V0 are visited, and once that happens, vertices in Vd+1 may start to be visited. This pattern repeats for all the remaining sets, that is, vertices in Vq may only be visited after all vertices in Vp were visited, for all p, q ∈ P , q > p+d. The rule is illustrated in the example below. Suppose a problem has eight priority sets, namely V0, V1, ..., V7, and d is set to 1. The only sets that may be visited right after the origin are V0 and V1, and the d-relaxed 25 26 Formulations for the CTSP-d priority rule imposes that V2 may only be visited after all the vertices in V0 have been visited. It is worth emphasizing the cumulative aspect of the rule - vertices in set V3 may only be visited after all vertices in V0 and V1 have been totally visited, vertices in set V4 can only be accessed after all vertices in V0, V1 and V2 have been totally visited, and so on. There are several formulations for the TSP based on the combination of an as- signment problem with binary variables and a set of subtour elimination constraints (Öncan et al., 2009; Bektaş and Gouveia, 2014), which will be the basis for the formu- lations presented in this chapter. Since all the formulations are based on a common set of variable, constraints and parameters, all these common elements are defined first before presenting the formulations. • Parameter: cij: Cost associated with the travel from vertex i to vertex j. • Variable: xij = 1 when arc (i, j) ∈ A is included in the tour, or 0 otherwise. Considering that we will present different formulations for the CTSP with d-relaxed priority rule, a general formulation for the TSP is given considering a generic subtour elimination constraint, as well as a generic d-relaxed priority rule constraint: min ∑ (i,j)∈A cijxij (3.1) s. t. ∑ i∈V :(i,j)∈A xij = 1, ∀j ∈ V (3.2) ∑ j∈V :(i,j)∈A xij = 1, ∀i ∈ V (3.3) xij ∈ {0, 1}, ∀(i, j) ∈ A (3.4) solution does not contain subtours (3.5) solution satisfies the d-relaxed priority rule (3.6) In this formulation, the objective function (3.1) represents the goal of minimizing trav- eling costs in the problem, given by the sum of the costs of the arcs selected for the tour. Constraints (3.2) and (3.3) define the flow in the solution, imposing that every vertex must be visited and left exactly once and constraints (3.4) define the domain of the variables xij. The other two sets of constraints, namely (3.5) and (3.6), are generic constraints representing the final aspects that a formulation for the CTSP-d must satisfy. Constraints (3.5) represent the subtour elimination constraint set, and so expressions (3.2)-(3.5) define a general model for the TSP, while constraints (3.6) are the d-relaxed priority rule constraints, which turns the classic TSP into the CTSP-d. Miller-Tucker-Zemlin based formulations 27 3.1 Miller-Tucker-Zemlin based formulations A simple way to obtain a compact formulation for the CTSP-d is adopting con- straints based on the classic Miller-Tucker-Zemlin (MTZ) subtour elimination con- straints (Miller et al., 1960). An important feature of this approach is it having a small number of variables and constraints, being an O(n2) formulation for variables and constraints. It is the basis for the formulations proposed in Panchamgam (2011) and Hà et al. (2020) for the CTSP-d, and also in Doan et al. (2021) for the VRP with d-relaxed priority rule. In its original form, the MTZ constraints are as in (3.7) ui − uj + nxij ⩽ n− 1, ∀(i, j) ∈ A, j ̸= 0, (3.7) where uj ∈ R+, j ∈ V is a set of continuous variables that represents the order of each vertex j in the tour. Despite the fact that uj are defined as real variables, constraints (3.7) ensure that their values will be integer values for any feasible solution. Since the variables uj represent the order of the vertices in the solutions, they can be used to define a d-relaxed priority rule constraint by doing ui + 1 ⩽ uj, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. (3.8) • Hence, the first formulation of the CTSP-d, denoted by MTZ1, is defined by (3.1)-(3.4) together with constraints (3.7), that guarantees subtour elimination, and the d-relaxed priority rule constraints (3.8). A way to strengthen this formulation is to add valid inequalities. In Hà et al. (2020), the authors present four sets of constraints, some of them proposed by Panchamgam (2011), that can be used for this purpose, which are given in (3.9)-(3.12). ∑ i∈Vp ∑ j∈Vq : (i,j)∈A xji = 0, ∀p, q ∈ P, q > p+ d (3.9) ∑ j∈Vp x0j = 0, ∀p ∈ P, p > d (3.10) ∑ i∈Vp xi0 = 0, ∀p ∈ P, p < Pmax − d− 1 (3.11) ∑ i∈Vp ∑ j∈Vq : (i,j)∈A xij ⩽ 1, ∀p, q ∈ P, q > p+ d (3.12) Since the vertices in a set Vq only can be visited after all the vertices in Vp have 28 Formulations for the CTSP-d already been visited for all q > p+ d, there will be no trips from vertices in Vq back to vertices in Vp, which is expressed in constraints (3.9). Moreover, a similar relationship can be observed for the origin vertex, which is not included in any priority set. This is ensured by constraints (3.10) and (3.11). The constraints in (3.12) do complement this relationship in the following way: when a vertex in Vp is left, the next vertex visited can belong to a Vq set such that q > p + d, but this can only occur once at most. It is worth emphasizing the constraints (3.9)-(3.11) are used to fix the values of sets of variables to zero, and this can be done in different ways. In this thesis we will consider the constraints as described in the text to maintain the same approach described in Panchamgam (2011) and Hà et al. (2020). Formulation MTZ1 was described in Hà et al. (2020) as well as another weaker formulation for the CTSP-d based on the Miller et al. (1960) constraint sets and the authors showed that adding (3.9)-(3.12) to MTZ-based formulations is an effective way to improve the formulation of the CTSP-d. • The formulation given by (3.1)-(3.4), (3.7)-(3.8) and (3.9)-(3.12) will be denoted by H2020. 3.1.1 Proposed valid inequalities for the MTZ1 formulation In this section, new constraint sets are proposed that may be included or reformu- lated to strengthen the MTZ1 formulation. Variable bounds and strengthening constraints (3.7) The first improvement proposed is to change the original MTZ formulation to its lifted version, proposed by Desrochers and Laporte (1991), where the constraints set (3.7) is reformulated as (3.13). ui − uj + nxij + (n− 2)xji ⩽ n− 1, ∀(i, j) ∈ A : j ̸= 0. (3.13) In Proposition 1, the constant bounds used in inequality (3.13) are analyzed and shown that they can be strengthened using the variable bounds the CTSP-d induce over uj. The first step to prove this result is to deduce these new bounds, which are themselves valid inequalities to the problem too. The d-relaxed priority rule imposes that for a fixed value p, all the vertices in priority set Vp must be visited before the vertices in Vq start to be visited, for all q > p + d. In this sense, the priority rule defines upper and lower bounds for the position that each vertex may appear in a sequence, depending on the value of parameter d: in the worst case, a vertex in Vp will figure in a solution Miller-Tucker-Zemlin based formulations 29 as soon as it is allowed by the rule and the last vertex in Vp will be right before the first vertex in Vq, q > p+ d. This condition can be represented using the cardinality of the priority sets, which defines the two following constraints sets: uj ⩾ p−d−1∑ r=0 |Vr|+1, ∀j ∈ Vp, d+ 1 ⩽ p ⩽ Pmax (3.14) uj ⩽ p+d∑ r=0 |Vr|, ∀j ∈ Vp, 0 ⩽ p ⩽ Pmax − d− 1 (3.15) Constraints (3.14) and (3.15) were not defined for all possible values of p since they would reduce to the original bounds of the variables, that is, they would turn into ui ⩾ 1 or ui ⩽ n− 1 in the non-included cases. But they can still be used to combine constraints (3.13), (3.14) and (3.15) to deduce tighter bounds to (3.13), as described in Proposition 1. Proposition 1. Constraints ui − uj + (Mij + 1)xij + (Mij − 1)xji ⩽ Mij, ∀(i, j) ∈ A, i, j ̸= 0 (3.16) are valid for the CTSP-d, where Mji = Mij = min{Pmax,q+d}∑ r=max{0,p−d} |Vr|−1, ∀i ∈ Vp, j ∈ Vq, p = 0, ..., Pmax, q = p, ..., Pmax, (3.17) and these are the best bounds for this constraints set. Proof. Inequalities (3.16) are a generalization of the lifted MTZ inequalities de- scribed in (3.13) and, for any feasible solution, the right-hand side of (3.16) is an upper bound for the difference ui − uj. In this sense, we need to prove the values in (3.17) represent upper bounds for this difference. First, notice that if Mij is an upper bound to ui−uj > 0, so it is an upper bound to uj −ui too since the last difference will be a negative value. It implies that Mij = Mji, ∀(i, j) ∈ A, i, j ̸= 0. Also, the bounds described in (3.14) and (3.15) can be extended to all the Vp sets in the following way: assume that, when p−d−1 < 0, the summation does not represent anything feasible so it can be ignored and expression (3.14) may be written as: uj ⩾ p−d−1∑ r=0 |Vr|+1, ∀j ∈ Vp, 0 ⩽ p ⩽ Pmax. Moreover, all the uj variables must satisfy the initial TSP upper bound given by uj ⩽ n− 1. In inequality (3.15), this is not included, since this is already a constraint 30 Formulations for the CTSP-d in the MTZ formulations. A simple way to denote the general case is: uj ⩽ min{p+d,Pmax}∑ r=0 |Vr|, ∀j ∈ Vp, 0 ⩽ p ⩽ Pmax. Now, let i ∈ Vp and j ∈ Vq such that q ⩾ p. The deduced bounds imply uj − ui ⩽ min{Pmax,q+d}∑ r=0 |Vr|− ( p−d−1∑ r=0 |Vr|+1 ) = min{Pmax,q+d}∑ r=max{0,p−d} |Vr|−1, ∀i ∈ Vp, j ∈ Vq, p = 0, ..., Pmax, q = p, ..., Pmax, (3.18) where the right hand side of (3.18) can be obtained by noticing that when p−d−1 < 0, the bound for ui is just 1, and so: uj − ui ⩽ min{Pmax,q+d}∑ r=0 |Vr|−1, ∀i ∈ Vp, j ∈ Vq, p = 0, ..., Pmax, q = p, ..., Pmax. (I) On the other hand, when p− d− 1 ⩾ 0, the terms Vr common to both the bounds of uj and ui will become 0 in the difference, and the result will be: uj − ui ⩽ min{Pmax,q+d}∑ r=p−d |Vr|−1,∀i ∈ Vp, j ∈ Vq, p = 0, ..., Pmax, q = p, ..., Pmax. (II) Combining the conclusions in (I) and (II), the bound is given by (3.18). It implies that the lowest value Mji may assume is the right hand side of the inequality. Extending this to Mij: Mji = Mij = min{Pmax,q+d}∑ r=max{0,p−d} |Vr|−1, ∀i ∈ Vp, j ∈ Vq, p = 0, ..., Pmax, q = p, ..., Pmax. (3.19) Noticing this relationship is valid for any i, j ∈ V , i ̸= j, i, j > 0, it gives a new stronger bound to the lifted MTZ inequalities for the CTSP-d. To see that this is the best possible bound for this constraints set, suppose that a smaller bound Bji ∈ R exists satisfying the condition Bji < Mji for some arc (j, i) ∈ A, and so, without loss of generality, uj ∈ Vq and ui ∈ Vp such that uj > ui. The definition of the problem allows a solution to be built such that: uj = q+d∑ r=0 |Vr| and ui = p−d−1∑ r=0 |Vr|+1. Miller-Tucker-Zemlin based formulations 31 In this case, the solution obtained satisfies uj − ui = Mji. But we are supposing the difference uj − ui ⩽ Bij < Mji, which is a contradiction. It proves that the smallest values the parameters Mij may assume in inequality (3.16) are the ones given in (3.19), which concludes the proof. Beyond the constraints sets proposed in (3.14) and (3.15), there are other known inequalities that may be used to deal with uj variables bounds. In the paper by Desrochers and Laporte (1991) there are three other constraints sets presented to lift the lower and upper bounds of the uj variables, namely: ui ⩾ 1 + ∑ (j,i)∈A:j ̸=0 xji, ∀i ∈ V, i ̸= 0 (3.20) ui ⩽ n− (n− 2)x0i − ∑ (i,j)∈A:j ̸=0 xij, ∀i ∈ V : (0, i) ∈ A (3.21) ui ⩾ (n− 2)xi0 + ∑ (j,i)∈A:j ̸=0 xji, ∀i ∈ V : (i, 0) ∈ A (3.22) The constraints in (3.20) stand for the fact that ui ⩾ 1, ∀i ∈ V , i ̸= 0, which is valid for both TSP and CTSP-d. In this case, however, constraints (3.14) grant greater bounds for uj such that j ∈ Vp, p ⩾ d + 1, and so constraints (3.20) can be restricted to the first d sets, which leads to the constraint set (3.23). ui ⩾ 1 + ∑ (j,i)∈A:j ̸=0 xji, ∀i ∈ Vp, 0 ⩽ p ⩽ d. (3.23) Constraints (3.21) and (3.22) refer to the vertices that come right before or after the origin in the tour. Considering the constraint sets (3.9)-(3.11) are part of the new model, these constraints can be added to the CTSP-d formulation in its original for- mulation, since values such that xi0 = 0 or x0i = 0 are already fixed to be zero. Finally, the value of u0 can be fixed in the following way: u0 = 0. (3.24) Strengthening constraints (3.8) In a similar way, the xij can be added to constraints set (3.8). In fact, the condition of the d-relaxed priority rule imposes that vertices in Vp sets with lower values of p must be visited before the vertices in a Vq set, where q > p+ d and, in the worst case, ui−uj = −1, for some i ∈ Vp, j ∈ Vq where q > p+d, since variables uk are integer for all feasible solutions. The fact that allows xij to be added to the constraint set is that if xij = 0, then uj − ui ⩾ 2, that is, vertex j will be at least 2 positions after vertex i 32 Formulations for the CTSP-d in the solution. Noticing that xji will never happen in (3.8), since vertices in Vq must succeed vertices in Vp, the new constraint set is ui + 2− xij ⩽ uj, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. (3.25) A more formal way to obtain this conclusion is by applying lifting to constraints (3.8), which are presented in Proposition 2. Proposition 2. Constraints (3.25) are valid for the CTSP-d. Proof. Constraints (3.25) are a lifted version of the original inequalities described in (3.8). The lifting process is applied to constraints (3.8) by adding new real parameters α, β and γ to the original constraints set, which leads to the following: ui + 1 + α + βxij + γxji ⩽ uj, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. The case where α = β = γ = 0 is the original set (3.8) of inequalities. We need to find the largest possible values for the parameters such that the new constraints set remains valid for the problem. As already discussed, we do not need to consider the γ parameter, since xji will always be zero for i ∈ Vp, j ∈ Vq when q > p + d, and so the final problem is to find the largest values for α and β in expression (3.26) ui + 1 + α + βxij ⩽ uj, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d, (3.26) In the case xij = 0: ui + 1 + α ⩽ uj, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. But uj > ui by hypothesis, and moreover uj −ui > 1, since xij ̸= 1, which implies that the largest value α can assume is α = 1. Considering the case xij = 1: ui +1+α+ β ⩽ uj ⇒ 1+α+ β ⩽ 1⇒ β ⩽ −α, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. Applying the value α = 1, the greatest value β can assume is β = −1. Replacing the values α = 1 and β = −1 in (3.26) gives the constraints set (3.25). Other inequalities Finally, two other valid inequalities sets are proposed to be included in the new MTZ-based formulation that are not directly related to the d-relaxed priority rule, but to general TSP properties. The first one follows from the fact the variables uj, j ̸= 0, Precedence variable-based formulations 33 will always map onto a permutation of the values 1, 2, ..., n − 1, and so constraint (3.27) must always be satisfied n−1∑ j=0 uj = n−1∑ j=0 j = n(n− 1) 2 . (3.27) The other relationship is that each arc on the graph can only be traversed in one direction in any feasible solution, and so the following two-vertices clique constraint set can be added: xij + xji ⩽ 1, ∀(i, j) ∈ A. (3.28) • The formulation MTZ2 assembles constraint sets (3.1)-(3.4), (3.9)-(3.12), (3.14)- (3.15), (3.16) applying the bounds presented in (3.17) and expressions (3.21)- (3.25) and (3.27)-(3.28). It is worth noticing that, like the MTZ1, this is a model with O(n2) variables and constraint sets. 3.2 Precedence variable-based formulations Besides the models proposed for the CTSP-d problem based on the classic MTZ formulation, there are other ways to formulate the problem and, in this section, we propose three new formulations based on precedence variables, which are more general than the ordering variables used in the MTZ formulations and can be used to build stronger models. In the whole section, the models will be formulated using yij variables as described below: yij = 1 if vertex i precedes (not necessarily immediately) vertex j in the tour, and 0 otherwise. As discussed in Sarin et al. (2005), the yij variables defined in this way allow precedence constraints just by setting a value of one to the variables matching the precedence condition of the problem, that is, adding the following set of constraints to the model: yij = 1, ∀j ∈ V \{0}, ∀i ∈ SPCj, (3.29) where SPCj is the set of vertices that must precede j in the tour. In the case of the CTSP-d problem, where the precedence relationships are defined through the Vp priority sets, the equations set (3.29) can be formulated in the following way: yij = 1, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. (3.30) The formulations proposed in this section are based on the formulations presented in the 34 Formulations for the CTSP-d review by Öncan et al. (2009) and in the papers in which each one of the formulations was originally proposed. All these formulations are proven to be stronger than the MTZ formulation for the classic TSP. Moreover, the proposed formulations may have valid inequalities included or some other kinds of reformulation techniques, which will be described here. Finally, we highlight that there are many other formulations for the TSP, and it is not our intention to exhaust the subject, but to demonstrate different ways to formu- late the CTSP-d, and present how the problem can be identified with the Precedence Constrained Traveling Salesman Problem (PCTSP). In addition to the formulations we adapted, it is an important research subject to study different ways to model the problem and how different formulations can be used to reach the best results in com- putational experiments. Moreover, it helps to recognize the CTSP-d as a particular or general case of other known variants of the TSP, which can lead to efficient solu- tion approaches to the problem. In this context, some recent papers which may be interesting are Godinho et al. (2014), where the authors present formulations for the Time-Dependent Traveling Salesman Problem, for which the asymmetrical TSP is a particular case, and then show how to develop an stronger formulation using precedence variables, Balma et al. (2018), where the authors work with strong multi-commodity flow formulations for the classical asymmetrical TSP and for the precedence-constrained case, and Gouveia et al. (2018), in which combinations and projections of flow models are used to improve the PCTSP. Formulation based on Gouveia and Pires (1999) The first formulation that is adapted to model the CTSP-d was based on precedence constraints proposed by Gouveia and Pires (1999), where the relationship between the variables xij and yij are defined by the constraints: xij − yij ⩽ 0, ∀i, j ∈ V : i, j > 0, i ̸= j (3.31) xij + yji ⩽ 1, ∀i, j ∈ V : i, j > 0, i ̸= j (3.32) xji + xij + yki − ykj ⩽ 1, ∀i, j, k ∈ V : i, j, k ̸= 0, i ̸= j ̸= k (3.33) xkj + xik + xij + yki − ykj ⩽ 1, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k (3.34) Constraints (3.31) represent the fact that if j is the successor of i in the tour, then i precedes j in the tour, and conversely, if i does not precede j in the tour, then j cannot be the successor of i. Similarly, constraints (3.32) ensure that if j is the successor of i, then j does not precede i, as well as if j precedes i, so we know that j is not the successor of i in the tour. Precedence variable-based formulations 35 Constraints (3.33) and (3.34) are lifted versions of the constraints: xij + yki − ykj ⩽ 1, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k, (3.35) which ensure the following condition: when the arc (i, j) is included in the tour, that is, when j is the successor of i in the tour, then yki ⩽ ykj, for all possible k ̸= 0, which means that every vertex k that precedes i in the tour must necessarily precede j. The purpose of using (3.33) and (3.34) instead of (3.35) is to build a stronger formulation (Öncan et al., 2009). One disadvantage of doing this is the fact that the formulation is bigger. It is important to notice that the combination of the expressions (3.1)- (3.4), (3.30), (3.31)-(3.32) and one of the constraints (3.33), (3.34) or (3.35) defines a formulation for the CTSP-d. • In this work we will call GP1 the formulation defined by (3.1)-(3.4), (3.30), and (3.31)-(3.34), which has O(n2) variables and O(n3) constraints. Formulation based in Sarin et al. (2005) The next formulation to model the CTSP-d is based on the formulation proposed by Sarin et al. (2005), where the variables xij and yij are related by the following constraints: yij ⩾ xij, ∀i, j ∈ V : i, j > 0, i ̸= j, (3.36) yij + yji = 1, ∀i, j ∈ V : i, j > 0, i ̸= j (3.37) yij + xji + yjk + yki ⩽ 2, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k, (3.38) x0j + xj0 ⩽ 1, ∀j ∈ V, j > 0 (3.39) The constraints set (3.36) ensures that if a vertex j is the successor of vertex i in a tour, then the associated variable yij is equal to 1, and at the same time it imposes that if j does not succeed i in the tour, then xij = 0. Another basic relationship between a pair of distinct vertices i and j in the graph is that j must be before or after i, condition that is guaranteed by constraints set (3.37). The set in (3.38) is a lifted version of yij + yjk + yki ⩽ 2, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k (3.40) and, together with (3.36) and (3.37), is sufficient to completely define the yij variables main condition, that is, it ensures yij = 1 if, and only if, j is successor of i in the tour. The constraints set (3.39) is not necessary to define the formulation, but they are included as valid inequalities to strengthen the formulation. 36 Formulations for the CTSP-d • In summary, an alternative formulation for the CTSP-d, which is denoted by SSB1, is defined by (3.1)-(3.4), (3.30) and (3.36)-(3.39), and it hasO(n2) variables and O(n3) constraints. Formulation based in Sherali et al. (2006) The last formulation adapted to fit the d-relaxed priority rule was proposed by Sherali et al. (2006). The authors deduce a reformulation for a model defined by (3.1)- (3.4), (3.30), (3.36)-(3.37) and (3.41)-(3.43), where the last constraints are defined as follows: yij ⩾ x0i, ∀i, j ∈ V : i, j > 0, i ̸= j, (3.41) yji ⩾ xi0, ∀i, j ∈ V : i, j > 0, i ̸= j, (3.42) −(1− xik) ⩽ yij − ykj ⩽ 1− xik, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k, (3.43) where inequalities (3.41) and (3.42) are a disaggregated version of (3.39) (notice that summing up the two expressions imply the same constraints as in (3.39) as a conse- quence of (3.37)), and inequalities (3.43) ensure that when xik = 1 for some i, k ∈ V , i, k > 0, i ̸= k, than the relation yij = ykj will be valid, for all j ∈ V , j > 0, i ̸= j ̸= k. In other words, the successors of i different from k will be successors of k too, and the same applies to the predecessors of k different from i. The new formulation is obtained by applying a Reformulation-linearization tech- nique (RLT) to the following non-linear relations: (0 ⩽ ykj ⩽ 1) · xik, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k,[ n−1∑ k=0,k ̸=i xik = 1 ] · yij, ∀i, j ∈ V : i, j > 0, i ̸= j, [ n−1∑ i=0,i ̸=k xik = 1 ] · ykj, ∀j, k ∈ V : j, k > 0, j ̸= k. To deduce the new constraints sets for the problem based on the non-linear constraints described above, a set of auxiliary variables must be defined representing the product of the variables xik and ykj in the following way: tkij = xikykj, ∀i, j, k ∈ V : i, j, k > 0, i ̸= j ̸= k. (3.44) Precedence variable-based formulations 37 After applying the RLT, the new constraints obtained are: 0 ⩽ tkij ⩽ xik, ∀i, j, k ∈ V : i, j, k ̸= 0, i ̸= j ̸= k (3.45)∑ k∈V \{0}: k ̸=i,j tkij + xij = yij, ∀i, j ∈ V : i, j > 0, i ̸= j (3.46) x0k + ∑ i∈V \{0}: i ̸=j,k tkij = ykj, ∀j, k ∈ V : j, k > 0, j ̸= k (3.47) The details involved in obtaining these new linear constraints can be found in the paper by Sherali et al. (2006), and all the steps are not reproduced here since this would not contribute to the goal of relating the d-relaxed priority rule and the already known formulations for the TSP. • The formulation is defined by (3.1)-(3.4), (3.30), (3.37)-(3.38), (3.41)-(3.42), (3.45)-(3.47) and is denoted by SST1. Note that this formulation is bigger than all the formulations previously described, since it has O(n3) variables and O(n3) constraints. 3.2.1 Proposed valid inequalities for the precedence variables- based formulations In both ways of formulating the TSP, i.e., using the position variables uj or the precedence variables yij, our main goal is to have variables representing global rela- tionships among the vertices in a tour, which allows the avoiding of subtours and the formulation of more complex variants of the problem, for example the CTSP-d dis- cussed in this thesis. How variables uj and yij are related is presented next, showing how the valid inequalities proposed for the formulation MTZ1 in Section 3.1.1 can be adapted and used to strengthen the GP1, SSB1 and SST1 formulations. The first point to emphasize is that any feasible solution is a tour where the origin vertex always figures before and after all the other vertices. The MTZ formulation follows the convention u0 = 0, which is ensured by the subtour elimination constraints. In the precedence-based formulations, the constraints do not define values for the yij variables, and so two new sets of constraints are imposed as follows: yi0 = 0, ∀i ∈ V \{0} (3.48) y0j = 1, ∀j ∈ V \{0} (3.49) These two sets of fixed values for yij allow an analogy between the way variables uj 38 Formulations for the CTSP-d relate to the feasible solutions in the MTZ-based formulations and the way variables yij relate to them in the precedence-based formulations. Constraints set (3.49) represents the fact that the vertex 0 is always the first one in any tour and it plays the same role as constraint (3.24) for the uj variables. Since the variables uj cannot represent the fact that 0 is the final vertex in the tour too (as long as u0 cannot assume two different values in the same solution), (3.48) is added to the precedence variable models and so the following relationship is valid: uj = ∑ k∈V : (k,j)∈A ykj, ∀j ∈ V, j ̸= 0. (3.50) In other words, the values variables uj assume, for each j ∈ V , represent how many vertices there are before j in a solution (except for the final vertex 0). Having this simple relationship well defined, any valid inequality proposed for the MTZ1 formulation can be directly adapted to any of the GP1, SSB1 or SST1 formulations through a direct substitution. Analyzing expressions (3.9)-(3.12) Another fact that can be used to strengthen the precedence variable models is that constraints sets (3.9)-(3.12) only rely on variables xij, which are used in all the formulations presented, and so they can be directly added to GP1, SSB1 and SST1. Proposition 3 shows that equation (3.9) does not need to be added into SSB1 and SST1. Proposition 3. Constraints (3.9) are redundant in SSB1 and SST1. Proof. Let i and j be two vertices such that i ∈ Vp and j ∈ Vq, where q > p + d. Constraints (3.30) imply yij = 1, and so yji = 0 by (3.37). But this conclusion leads to xji = 0, as a consequence of (3.36). And since i and j are arbitrary, xji = 0, for any i ∈ Vp, j ∈ Vq such that q > p+ d. Summing all these variables gives ∑ i∈Vp ∑ j∈Vq : (i,j)∈A xji = 0, ∀p, q ∈ P, q > p+ d, which is exactly the result in (3.9), and so these constraints would not make any difference to models SSB1 and SST1 if added. For the GP1 formulation, the condition in (3.37) is not ensured by its constraints sets, and so the conclusion reached for the two other precedence variables formulations is not valid. But since this condition must be valid for any feasible solution, add (3.37) Precedence variable-based formulations 39 to the new CTSP-d formulation based on GP1 and the result follows, since constraint sets (3.36) and (3.31) are the same. Moreover, in a formulation including (3.31) and (3.37), (3.32) is not needed, since it will follow from a direct substitution of (3.37) in (3.31). Analyzing variables bounds (3.14)-(3.15) and (3.21)-(3.23) Constraint sets (3.14) and (3.15) can be adapted to produce lower and upper bounds to the sum of the yij variables in the following way: ∑ k∈V : (k,j)∈A ykj ⩾ p−d−1∑ r=0 |Vr|+1, ∀j ∈ Vp, d+ 1 ⩽ p ⩽ Pmax (3.51) ∑ k∈V : (k,j)∈A ykj ⩽ p+d∑ r=0 |Vr|, ∀j ∈ Vp, 0 ⩽ p ⩽ Pmax − d− 1 (3.52) However, this would not lead to a better formulation, which is proved in Proposition 4. Proposition 4. Constraints (3.51) and (3.52) are redundant in the precedence-variable models. Proof. In fact, for each fixed j ∈ V \{0}, j ∈ Vp, for some p value, and so it follows from (3.30), (3.37) and (3.49) that ∑ k∈V : (k,j)∈A ykj = y0j + ∑ k∈Vr: r=0,...,p−d−1 ykj + ∑ k∈Vr: r=p−d,...,p+d ykj + ∑ k∈Vr: r=p+d+1,...,Pmax ykj = = 1 + p−d−1∑ r=0 |Vr|+ ∑ k∈Vr: r=p−d,...,p+d ykj ⇒ 1 + p−d−1∑ r=0 |Vr|⩽ ∑ k∈V : (k,j)∈A ykj ⩽ p+d∑ r=0 |Vr|. The lower and upper bounds found are exactly the same as in (3.51) and (3.52), and so the constraints are redundant for the models, as stated. Identity (3.50) can be applied to (3.21), (3.22) and (3.23), which leads to new constraints sets: ∑ k∈V : (k,i)∈A yki ⩽ n− (n− 2)x0i − ∑ (i,j)∈A:j ̸=0 xij, ∀i ∈ V : (0, i) ∈ A (3.53) 40 Formulations for the CTSP-d ∑ k∈V : (k,i)∈A yki ⩾ (n− 2)xi0 + ∑ (j,i)∈A:j ̸=0 xji, ∀i ∈ V : (i, 0) ∈ A (3.54) ∑ k∈V : (k,i)∈A yki ⩾ 1 + ∑ (j,i)∈A:j ̸=0 xji, ∀i ∈ Vp, 0 ⩽ p ⩽ d. (3.55) Adapting expressions (3.27) and (3.28) Note that (3.28) is a generalization of the result in (3.39) for all the vertices i ∈ V and can be added to the models, since it only depends on the xij variables. But this would not produce any improvement in SSB1 and SST1, as shown in Proposition 5. Proposition 5. Constraints (3.28) are redundant for SSB1 and SST1. Proof. Both the SSB1 and SST1 formulations consider the constraints set (3.36). This implies that, for any i, j ∈ V , i ̸= j, i, j > 0, the following condition must hold:{ xij ⩽ yij xji ⩽ yji ⇒ xij + xji ⩽ yij + yji = 1, ∀i, j ∈ V : i, j > 0, i ̸= j (3.56) In formulation SSB1, the case where i = 0 or j = 0 is directly considered as the constraints set (3.39). Finally, for the formulation SST1, summing the constraints in (3.41) and (3.42) gives a similar conclusion to that in (3.56) and so the result is valid. Constraint set (3.37) is included in all precedence-variable reformulations and so the result in (3.56) will be valid for the GP1-based reformulation too. In this case, the only constraints that can be added to this new formulation are (3.39) or (3.41) and (3.42), which include the initial vertex. The last constraint set proposed for MTZ2 to be analyzed is (3.27). In this case, applying the (3.50) relationship derives valid inequalities to the precedence-variable formulations. This gives the result presented in (3.57). ∑ (i,j)∈A yij = n(n− 1) 2 . (3.57) Since subtour elimination is already part of the original models, and the d-relax priority rule is defined in the models through (3.30), a summary of the new formulations for the precedence-variable formulations can be presented. • (GP2) is given by (3.1)-(3.4), (3.10)-(3.12), (3.30), and (3.31), (3.33)-(3.34), (3.37), (3.39), (3.53)-(3.55), (3.57). Precedence variable-based formulations 41 • (SSB2) is defined by (3.1)-(3.4), (3.10)-(3.12), (3.30) and (3.36)-(3.39), (3.53)- (3.55), (3.57). • (SST2) is given by (3.1)-(3.4), (3.10)-(3.12), (3.30), and (3.37)-(3.38), (3.41)- (3.42), (3.45)-(3.47), (3.53)-(3.55), (3.57). The size order of the formulations, the section in which each formulation is de- scribed, as well as its reference numbers are shown in Table 3.1. Besides the fact that adding valid inequalities increases the formulation size, the size order of the formula- tions is still the same for all of them. 42 Formulations for the CTSP-d V ariables C onstraints Form ulation H 2020 O (n 2) O (n 2) (3.1)-(3.4),(3.7)-(3.12) M T Z1 O (n 2) O (n 2) (3.1)-(3.4),(3.7)-(3.8) M T Z2 O (n 2) O (n 2) (3.1)-(3.4),(3.9)-(3.12),(3.14)-(3.17),(3.21)-(3.25),(3.27)-(3.28) G P 1 O (n 2) O (n 2) (3.1)-(3.4),(3.30)-(3.34) G P 2 O (n 2) O (n 2) (3.1)-(3.4),(3.10)-(3.12),(3.30)-(3.31),(3.33)-(3.34),(3.37),(3.39),(3.53)-(3.55),(3.57) SSB 1 O (n 2) O (n 3) (3.1)-(3.4),(3.30),(3.36)-(3.39) SSB 2 O (n 2) O (n 3) (3.1)-(3.4),(3.10)-(3.12),(3.30),(3.36)-(3.39),(3.53)-(3.55),(3.57) SST 1 O (n 3) O (n 3) (3.1)-(3.4),(3.30),(3.37)-(3.38),(3.41)-(3.42),(3.45)-(3.47) SST 2 O (n 3) O (n 3) (3.1)-(3.4),(3.10)-(3.12),(3.30),(3.37)-(3.38),(3.41)-(3.42),(3.45)-(3.47),(3.53)-(3.55),(3.57) Table 3.1: Size order and reference num bers ofeach form ulation 4 Formulations and a heuristic approach for the CluVRP-d The aim of the VRP-RPR, considered in this thesis, is to find up to Gmax vehicle routes starting and ending at the depot, such that: each vertex is visited exactly once and its demand (dj) is satisfied; the total demand quantity of any route does not exceed the vehicle capacity capk; the d-relaxed priority rule is satisfied (locally or globally); a weighted summation of the total travel cost (cijk) and fixed cost (Ck) is minimized. In this chapter we will call the VRP-PRP version considered here as the Clustered Vehicle Routing Problem with d-relaxed priority rule (CluVRP-d). The extension of the d-relaxed priority rule in the context of VRP can be divided in two approaches, namely the Local Timing model and the Global Timing model. The first case refers to the application of the rule on the tours described by each working vehicle in isolation, and so the rule is obeyed by each vehicle only over the areas that vehicle visits. However, this rule does not consider any global hierarchy constraints; that is, once a vehicle finishes its work on the higher priority areas in which it was assigned to work, the vehicle may go to the lower priority areas in its tour, which may occur regardless of the existence of high priority areas still not visited in the other vehicles routes. On the other hand, the Global Timing model requires that the d-relaxed priority rule must be obeyed not only by each vehicle individually, but it includes a global hierarchy constraint across all the tours; that is, a vehicle can only visit a lower priority area after all the higher priority areas have been visited by some vehicle. That is, the priority rule must be observed across the fleet of vehicles. Both the Local and the Global Timing Models were initially described in Pan- chamgam (2011), and the Local Timing case was treated in Doan et al. (2021) too. In both cases the main aspects described are: each vertex may be visited, at most, once, each vehicle has a maximum capacity and there are maximum lengths for the routes of each vehicle. Considering these features, the priorities must satisfy two pri- ority rules, namely the “d-relaxed priority rule” and the “order of demand fulfillment (ODF)”, which states that “at the end of the planning horizon, all the higher priority vertices must be visited before the vehicles can move to any lower priority vertex”. The 43 44 Formulations and a heuristic approach for the CluVRP-d objective considered is to minimize the total traveling costs in the tours, as well as minimize the total demand not met in the vertices not visited. An important property proposed by Doan et al. (2021) is to present a problem where not all areas must be visited and, in this case, the d-relaxed priority rule allows tours that strongly weaken the original priorities hierarchy. It can be seen that this rule may allow the construction of tours in which a vehicle will only visit vertices with medium or lower priorities, mixing the different priorities requirements in a very aggressive approach. Moreover, supposing the different priority levels have relevant practical implications, to accept solutions focusing on lower priority vertices instead of the higher priority ones may lead to solutions that cannot be used in real world contexts. It motivates the definition of the ODF rule. According to Panchamgam (2011) and Doan et al. (2021), the ODF imposes that at the end of the planning horizon, if a vehicle has sufficient capacity, then all vertices with higher priorities must be visited before the vehicle be allowed to visit a vertex with lower priority level. In other words, combining the d-relaxed priority rule and the ODF in a Vehicle Routing Problem in which not all vertices must be visited ensures that, in the solution routes, at least the areas with higher priority levels were visited. We emphasize that it does not extinguish the possibility of a vehicle visiting only areas with medium or lower priorities. However, we can impose that the other vehicles visit the areas with higher priorities, not allowing the construction of routes that “ignore” the existence of priorities to reach smaller/cheaper routes or visits only to vertices with lower demands. In this study we propose the Clustered Vehicle Routing Problem with d-relaxed priority rule (CluVRP-d) as a direct generalization of the Clustered Traveling Salesman Problem with d-relaxed priority rule (CTSP-d), described in Panchamgam (2011), Hà et al. (2020), Doan et al. (2023) and Teixeira and de Araujo (2024), which have stricter constraints than the VRP-RPR studied in Doan et al. (2021). In particular, we impose that all vertices in the problem must be visited, there is a capacity for the vehicles but not a maximum route length, and the d-relaxed priority rule must be valid in any feasible solution. Requiring that all vertices must be visited removes any consequences of not visiting a vertex, a possibility which was considered by Doan et al. (2021). It is also sufficient to ensure that the solutions do not skip the priority levels, and so removes the need to consider the ODF rule. Formulations for the CluVRP-d 45 4.1 Formulations for the CluVRP-d Both the proposed models (Local Timing and Global Timing) are based on a com- mon vehicle routing problem structure, which we present below, leaving the details related to each approach for its own section. • Sets: G = {1, 2, ..., Gmax}: Set of available vehicles. • Parameters: dj: Total demand of vertex j ∈ C. Ck: Fixed cost of vehicle k ∈ G; capk: Total demand vehicle k ∈ G can meet over the planning horizon (total working capacity of vehicle k ∈ G); cijk: Cost to move vehicle k ∈ G from vertex i ∈ V to vertex j ∈ V ; • Variables: xijk = 1 when arc (i, j) ∈ A is included in the tour assigned to vehicle k ∈ G, or 0 otherwise. yk = 1 when there is work assigned for vehicle k ∈ G, or 0 otherwise. • Formulation: min ∑ (i,j)∈A ∑ k∈G cijkxijk + ∑ k∈G Ckyk (4.1) ∑ i∈V : (i,j)∈A ∑ k∈G xijk = 1, ∀j ∈ C (4.2) ∑ i∈V : (i,h)∈A xihk = ∑ j∈V : (h,j)∈A xhjk, ∀h ∈ C, k ∈ G (4.3) ∑ (i,j)∈A djxijk ⩽ capk, ∀k ∈ G (4.4) ∑ j∈V :(0,j)∈A x0jk = yk, ∀k ∈ G (4.5) The function defined in (4.1) represents the objective of minimizing the logistics costs involved in the problem, given by the sum of the traveling costs between the vertices and the fixed costs associated with the use of each vehicle in the fleet. The feasible solutions for the problem are given by a set of disjoint tours which start at vertex 0 and such that every vertex different from 0 is visited only once. The constraints in (4.2) ensure the visiting condition of the vertices, while the constraints in (4.3) define the flow condition on the problem - a vertex can be visited by a vehicle if, 46 Formulations and a heuristic approach for the CluVRP-d and only if, it is then left by this same vehicle. We consider in the formulations the application of the d-relaxed priority rule in the capacitated vehicle routing problem with heterogeneous fleet, that is, each vehicle k ∈ G has its own maximum working capacity capk associated, which must not be exceeded. The constraints in (4.4) ensure this condition. The variables yk represent the decision about using or not the available vehicles k ∈ G. The constraints in (4.5) impose that vehicle k ∈ G leaves the origin to visit some vertex if, and only if, vehicle k was chosen for the assignment of some tour, that is, the associated variable yk has a value of 1. 4.1.1 Local Timing Model The Local Timing Model refers to the case where the d-relaxed priority rule is applied to the route of each working vehicle in isolation, which means that the rule must only be satisfied by each vehicle over the vertices it was assigned to visit. In practical contexts, this can imply that once a vehicle finishes its work on the higher priority vertices in which it was assigned, it may move to the lower priority vertices it needs to visit, even if there are other high priority vertices not visited in the other vehicle tours. In other words, the vehicles must only satisfy the d-relaxed priority rule locally on their own routes, which motivates the name of the model. To propose a formulation for the Local Timing Model, we will add a set of auxiliary position variables to the common formulation (4.1)-(4.5) to model the d-relaxed priority rule. Here, we describe the new variables as well as the constraints relating them to the original xijk variables. • New variable set: ujk: Variable representing the position of vertex j in the tour of vehicle k when vertex j is assigned to the vehicle k tour, or 0 otherwise. Since we are supposing not all the vehicles need to be used in a solution but they all are available in the depot, we may fix u0k = 0, for all k ∈ G. In this case, the first constraints we consider are 0 ⩽ ujk ⩽ n− 1, ∀j ∈ V, k ∈ G (4.6) u0k = 0, ∀k ∈ G (4.7) • Visiting order constraints/Sub-tour elimination: as described above, variables ujk represent the order of visiting each vertex in the tours of the working vehicles. Con- straints (4.8) are the lifted version of the classical Miller et al. (1960) constraints, which was proposed for the Traveling Salesman Problem by Desrochers and Laporte (1991), and they suffice to eliminate sub-tours in the problem. By adding the constraints Formulations for the CluVRP-d 47 set (4.9), we impose that when j immediately follows i in the tour of a vehicle, then ujk = uik + 1, that is, the variables represent the succession relationship. uik − ujk + nxijk + (n− 2)xjik ⩽ n− 1, ∀(i, j) ∈ A, k ∈ G : j ̸= 0 (4.8) uik − ujk − nxijk − (n+ 2)xjik ⩾ −n− 1, ∀(i, j) ∈ A, k ∈ G : i ̸= 0 (4.9) • Relationship between ujk and xijk: since we are considering that the variables ujk will represent the position of the vertex j in the tour of vehicle k, it can only be different from 0 when, in fact, the vertex j belongs to the tour of vehicle k, and so we need to consider the following constraints set: ujk ⩽ (n− 1) ∑ i∈V : (i,j)∈A xijk, ∀j ∈ C, k ∈ G. (4.10) • d-relaxed priority rule: our main goal is to apply the d-relaxed priority rule to the vehicle routing problem. In the case of the local timing model, it can be done as in equations (4.11), where we impose that when i, j ∈ C are included in the tour of a vehicle k, then they must satisfy the d-relaxed priority rule condition. uik + 1 ⩽ ujk + n 2− ∑ h∈V : (h,j)∈A xhjk − ∑ h∈V : (h,i)∈A xhik  , ∀i ∈ Vp, j ∈ Vq, k ∈ G : p, q ∈ P, q > p+ d. (4.11) 4.1.2 Global Timing Model Different from the Local Timing case, in the Global Timing Model the d-relaxed priority rule must be observed not only in the tours of each vehicle individually, but also across the tours; that is, a vehicle can only visit a lower priority vertex after all the higher priority vertices have been visited by some vehicle. That means all the vehicles must satisfy the d-relaxed priority rule globally in the problem, which leads to the name of the model. In the Global Timing Model, the auxiliary variables needed to formulate the prob- lem are cumulative time variables, which allows to measure when all the higher priority vertices were visited and allow the lower priority vertices to be visited. We describe below the extra variables and parameters used in the formulation, and the constraints proposed. • New parameters: 48 Formulations and a heuristic approach for the CluVRP-d tijk: Time vehicle k spends to travel between vertices i ∈ V and j ∈ V ; Mij: A set of big numbers used as auxiliary parameters to define the cumulative time variables; M : A sufficiently big number, representing the maximum value any cumulative time variable may assume. Since we are not imposing any total planning horizon timing for any vehicle, this value must be big enough to not prohibit any feasible tour. In the situation that the solutions should have an upper bound on the total time of the tours, this value can be used to impose it. • New variable: sjk: Variable representing the total cumulative traveling time of vehicle k from the depot node 0 up to vertex j, if j belongs to the tour of vehicle k, or 0 otherwise. • Cumulative time constraints/Sub-tour elimination: analogous to constraints (4.8) and (4.9) in the local timing model, here in the global timing model we have constraints based on the lifted version of the Miller et al. (1960) model to eliminate sub-tours and impose a relationship between a pair of variables sik and sjk when arc (i, j) ∈ A is included in the tour assigned to vehicle k. In this case, constraints (4.12) and (4.13) guarantee that when j succeeds i in a tour, the equality sjk ⩾ sik + tij must hold, and s0k = 0, for all k ∈ G. sik − sjk +Mijxijk ⩽ Mij − tij, ∀(i, j) ∈ A, k ∈ G : j ̸= 0 (4.12) s0k = 0, k ∈ G (4.13) • sik variables setup/tour time limit: the constraints presented in (4.14) ensure two important conditions to the model: when the vehicle k does not visit the vertex j, we will have sjk = 0, otherwise if the variables sjk have an upper bound, it can be imposed using the M parameter value. sjk ⩽ M ∑ i∈V : (i,j)∈A xijk, ∀j ∈ C, k ∈ G. (4.14) • d-relaxed priority rule: as discussed above, the difference between the Local Timing Model and the Global Timing Model is that, in the first case, the d-relaxed priority rule must be satisfied by each vehicle k ∈ G in an isolated manner, while in the global case the rule is obeyed by all the vehicles across the tours. In this sense, constraints (4.15) ensure the validity of the d-relaxed priority rule by imposing that given two vertices i, j ∈ V , if i ∈ Vp, j ∈ Vq and q > p+ d, then a vehicle can only visit A heuristic approach for the Local Timing case 49 j after the vertex i has already been visited. ∑ k∈G sik ⩽ ∑ k∈G sjk, ∀i ∈ Vp, j ∈ Vq : p, q ∈ P, q > p+ d. (4.15) 4.1.3 Valid Inequalities We included some common sets of valid inequalities in both the Local and Global Timing Models to build stronger formulations, which are presented below. ∑ i∈Vp ∑ j∈Vq : (i,j)∈A xjik = 0, ∀k ∈ G, p, q ∈ P : q > p+ d (4.16) ∑ i∈Vp ∑ j∈Vq : (i,j)∈A xijk ⩽ 1, ∀k ∈ G, p, q ∈ P : q > p+ d (4.17) ∑ k∈G (xijk + xjik) ⩽ 1, ∀(i, j) ∈ A. (4.18) Constraints in (4.16) represent the fact that the vertices in a set Vq can only be visited after all the vertices in the Vp sets have all been visited first, and it must hold for all q > p+ d. Moreover, for any vehicle k ∈ G, when a vertex belonging to a Vp set is left, it is possible that the next vertex in the tour belongs to a Vq set such that q > p+d, but it can only occur, at most, once, which is represented in (4.17). Both constraints (4.16) and (4.17) were proposed by Doan et al. (2021), as generalizations for the multi-vehicle case of the valid inequalities proposed by Hà et al. (2020) for the CTSP-d. Finally, the constraints set (4.18) represents the fact that each arc (i, j) can be included in the solution assigned to a single vehicle, and it can be traversed in only one direction, that is, if a vehicle k ∈ G goes from i to j, then it cannot go back from j to i. 4.2 A heuristic approach for the Local Timing case Considering the Local Timing case, we proposed a three-phase heuristic procedure based on clustering the vertices and after that, generating an initial feasible solution and then improving it through a variable neighborhood search method. This is a common approach for the Vehicle Routing Problem, as pointed out by Laporte et al. (2014). To cluster the vertices and produce the initial solution we solve a mixed integer linear programming problem based on the p-Medians problem, after that the variable neighborhood search procedure intra-route and inter-route moves are applied. The subject of metaheuristics for variants of the VRP is very comprehensive, including many different approaches based on local search methods, population based methods, 50 Formulations and a heuristic approach for the CluVRP-d matheuristics, as well as other kinds of solution procedures. A broader discussion can be found in Laporte et al. (2014) and Desaulniers et al. (2014). 4.2.1 Clustering Phase The goal of the clustering phase is to pack the vertices of the problem into clusters in a way that, for each cluster, it is possible to order its vertices and produce a valid initial route. It is not our intent to discuss the topic “clustering algorithms” but just show how to apply the method as part of the initialization phase of our heuristic. The reader can find deeper analysis about the algorithm in Ikotun et al. (2023). The clustering procedure is split in two steps that are described as follows. Uncapacitated clustering step The first step in the proposed heuristic is clustering the vertices according to their coordinates. To do so, we first use the K-means method (Lloyd, 1982), a classical and computationally efficient clustering algorithm, which can be applied in Python using the SciKit Learn (sklearn) module (https://scikit-learn.org/stable/) treating the data in the following way - the coordinates of the vertices are put in a table, where each vertex is in a line, the x coordinate is in the first column, while its y coordinate is in the second column. The K-means method is an iterative procedure based on clustering vertices using centroids that are updated at each iteration until the algorithm stops. The centroids of the clusters are initially created in determined positions and each vertex is assigned to the cluster with the nearest centroid position in a way that each vertex belongs to only one cluster. In the next step, for each cluster, the mean of its vertices features are calculated, the position of the centroid is changed to this value and after that the vertices are reassigned to the clusters with the nearest centroids again. In Algorithm 1 we describe the complete K-means method, summarizing its two main steps - assign the vertices to the clusters and move the centroids to their new position. The stopping criterion of the procedure is when it reaches a maximum number of iterations, but the method can be stopped when the clusters are no longer being improved after an iteration, or the change in the centroids positions are smaller than a given value. We highlight that there are different ways to initialize the clusters positions in the K-means method, and so the final result obtained by the method may vary accord- ingly. In our case we are initializing the clusters positions randomly, using the sklearn embedded functionalities, and so the results tend to be slightly different each time a problem is solved. The final result will be a set of K clusters that do not have any https://scikit-learn.org/stable/ A heuristic approach for the Local Timing case 51 Algorithm 1 K-Means procedure Input: A set of vertices {x1, ..., xn} ∈ Rm, a number of centroids K and a maximum number of iterations IterQty. Output: A set of disjoint clusters S1, ..., SK with centroids coordinates µ1, .., µK ∈ Rm. 1: t← 0 2: Randomly initialize µ0 1, ..., µ 0 K 3: repeat 4: t← t+ 1 5: for k = 1 to K do 6: St k ← ∅ 7: end for 8: for j = 1 to n do 9: k∗ ← argmin k=1,...,K {∥xj − µt k∥2} // Assign vertex xj to the nearest centroid 10: St k∗ ← St k∗ ∪ {xj} 11: end for 12: for k = 1 to K do 13: µt k ← 1 |St k| ∑ xj∈St k xj // Update centroids position 14: end for 15: until t ⩾ IterQty 16: for k = 1 to K do 17: Sk ← St k 18: µk ← µt k 19: end for vertices in common and all vertices are assigned to exactly one cluster. Moreover, since the positions of the clusters are defined iteratively using the mean of the coordinates of the vertices assigned to each of them, in the general case, the coordinates of the clusters do not need to correspond to the coordinates of the areas in the Routing Problem. It gives a clustering procedure for the problem vertices and, after this step, for each cluster, it is only necessary to assign a vehicle to its areas and the Routing Problem is reduced to a set of Traveling Salesman Problems. The impracticality of using this approach solely is that, by doing so, we are not taking into account the vehicles working capacity. By just clustering the vertices based on (x, y) coordinates some clusters might represent infeasible problems, including a number of vertices that would require a working capacity level that none of the vehicles available could meet. To overcome this, we propose a second step in the clustering procedure, which is described as follows. Capacitated clustering step The second step to cluster the vertices is to build a large number of clusters using only the coordinates of the vertices through the K-means method, as described pre- viously, and then to use the centroids coordinates as parameters in an optimization 52 Formulations and a heuristic approach for the CluVRP-d model based on the p-Medians Problem to merge these clusters, obtaining a smaller number of clusters needed to solve the Routing Problem, but imposing the working capacity as constraints in this model. The p-Medians problem is a classical integer programming problem in the field of discrete location that can be briefly described as follows: given a set of n users that must be supplied by exactly p facilities (usually called medians in the context of the problem), how to choose these facilities among a bigger set of m potential facilities in a way that minimizes the allocation costs of the users to the p medians. We can model this problem representing the m potential facilities and the n users as the vertices of a graph, and when a facility is chosen to be part of the solution, the vertices assigned to it are related to the facility by edges in the graph. A common case in the literature is to consider the set of users as the set of potential facilities, which is not the case we are considering. For more details about the problem the reader can refer to Marín and Pelegrín (2019). Our goal in this step of the heuristic is to choose a set of exactly r clusters repre- senting the r vehicles that will be used to solve the problem. Considering the vertices that must be visited in the Routing problem are the users in the p-Medians problem, by pre-clustering the vertices using the K-means method, we obtained a set of poten- tial facilities as the K centroids selected, and the present step is to select exactly r of these centroids aiming to minimize the sum of the distances among the vertices to the selected centroids and the fixed costs of the vehicles associated to them. It can be seen as an application of the p-Medians problem, where for the sake of notation we are representing the p parameter as r. To model the problem we will call the set of centroids selected by the K-means as centroids candidates and the centroids chosen to be part of the final solution as definitive centroids. To present our model based in the p-Medians problem, we need to define the fol- lowing auxiliary parameters and variables: • New set: K: Set of K centro