Efficient Matheuristics for Solving Production-Routing Problems

Computers Industrial Engineering 171 2022 108369 Available online 25 June 2022 03608352 2022 Elsevier Ltd All rights reserved Contents lists available at ScienceDirect Computers Industrial Engineering journal homepage wwwelseviercomlocatecaie Efficient matheuristics to solve a rich productionrouting problem Allexandre Fortes a Ricardo Camargo b Leandro Reis Muniz a Fátima Machado de Souza Lima c Fernanda dos Reis Cota d a Departamento de Engenharia Mecânica e Produção Universidade Federal de São João delRei Centro São João del Rei Minas Gerais Brazil b Departmento de Engenharia de Produção Universidade Federal de Minas Gerais Pampulha Belo Horizonte Minas Gerais Brazil c Departamento de Ciências Administrativas Universidade Federal de Minas Gerais Pampulha Belo Horizonte Minas Gerais Brazil d Departamento de Engenharia de Produção Administração e Economia Universidade Federal de Ouro Preto Bauxita Ouro Preto Minas Gerais Brazil A R T I C L E I N F O Keywords Productionrouting problem Iterated local search Hybrid methods Matheuristics A B S T R A C T We present a rich productionrouting problem having limited production and storage capacities at the plant limited storage capacity at the clients a heterogeneous fleet subjected to a maximum riding time and allowing for backorders to meet unfulfilled demands at penalty cost As the problem scales quickly with the number of customers periods products and vehicles three hybrid twolevel decomposition approaches using a topdown strategy were devised The top tier determines the production and inventory levels and the distribution of goods via CPLEX that is it makes tactical decisions while the bottom tier heuristically routes a heterogeneous fleet in each period that is it makes operational decisions The proposed methods rely on an iterated local search framework that combines tailored perturbation schemes prioritizing either tactical or operational decisions or both The main new feature of the algorithms is the adoption of an implicit cost that estimates the delivery routing costs when making production holding and transportation decisions This implicit cost serves as an important guide to obtain improved solutions The algorithms were tested over an extensive set of instances and the results demonstrated that all methods overcome CPLEX by obtaining more better and faster solutions with much less computational effort The devised heuristic which prioritizes operationallevel decisions during the perturbation phase attained the best overall results 1 Introduction The inventoryrouting problem IRP and its many variants combine inventory management and vehiclerouting decisions to minimize costs while improving supply chain performance and customer service levels When production decisions are also incorporated during modeling it gives rise to a more general problem known as the production routing problem PRP thereby enabling further savings The savings achieved by combining such decisions into a single problem have already been discussed by several authors addressing different practi cal applications eg Absi Archetti DauzèrePérès Feillet and Sper anza 2018 Archetti and Speranza 2016 Brown Keegan Vigus and Wood 2001 Çetinkaya Üster Easwaran and Keskin 2009 Côté Guastaroba and Speranza 2017 and Stålhane Andersson Chris tiansen and Fagerholt 2014 being Chandra 1993 one of the early works to point them out Most of these applications are within the context of vendormanaged inventory systems wherein a supplier has to efficiently manage retail ers stock levels by deciding the time and product quantities to deliver to the retailers so that stockouts do not occur Corresponding author Email addresses afortesufsjedubr A Fortes rcamargodepufmgbr R Camargo leandroreisufsjedubr LR Muniz fatimamslimaufmgbr FMS Lima fernandacotaalunaufopedubr FR Cota Stockout is a situation wherein an item is out of stock either by defective shelf replenishment practices or surges in the retailers demands that were not properly foreseen or expected by the supplier However as highlighted by Brahimi and Aouam 2016 a supply chain disruption due to limited production or depot storage capacities fleet related problems that cause delays in replenishment or as the whole world is now enduring the outbreak of a pandemic Furthermore there are many practical cases especially in the fast moving consumer goods FMCG industry where stockouts are dealt on a regular basis forcing managers to handle backorders that is managers have to satisfy a retailers demand for a given period on subsequent periods on a regular basis FMCG are goods having a useful life shorter than a year and are bought fairly frequently with recurring expenditures Examples of such products are food beverages personal care household and cleaning products apparel and footwear Despite its large market size the US FMCG market size alone was 635 billion in 2019 and despite the rich literature on PRP Adulyasak et al 2015b there are only a httpsdoiorg101016jcie2022108369 Received 10 November 2021 Received in revised form 2 May 2022 Accepted 21 June 2022 Computers Industrial Engineering 171 2022 108369 2 A Fortes et al few studies that have addressed backordering in the context of PRP namely Brahimi and Aouam 2016 and Bayley and Bookbinder 2019 Brahimi and Aouam 2016 proposed two mixedinteger linear pro grams for a rich PRP with a production facility with limited production and storage capacities considering multiple products delivered by a homogeneous fleet to geographically dispersed clients that also have limited storage capacities while allowing for backordering at a penalty cost They proposed a relaxandfix heuristic combined with a local search approach to solve their models and computationally assessed the impact and effects of different configurations on the production storage capacities as well as the demand patterns on backordering The results highlighted the importance of having backorders when addressing capacitated systems Finally Bayley and Bookbinder 2019 also proposed a branchandcut method for rich PRP with backorders but with no clear description of their algorithm and the generated instances Here we present a rich variant of PRP comprising deter mining at minimal total cost a production and distribution plan that fulfills periodic demands of multiple products by different clients while controlling the inventory levels at the production plant and clients and satisfying production and storage capacities that allows for back orders for the unfilled demands in a timely fashion and assumes a heterogeneous capacitated fleet with a maximum riding time to deliver the products in each period of the planning horizon The total cost is made up of production and distributionrelated costs The former includes the production setup holding and inventory costs while the latter accounts for vehicle fixed and routing costs To solve this new variant we designed three efficient hybrid it erated local search matheuristics that decompose the problem into two smaller subproblems one a tactical problem that decides what to make and to deliver at each period while satisfying the production and storage capacities and the other an operational problem that actually plans routes satisfying the maximum riding time of the fleet While the former subproblem is solved exactly with CPLEX the latter is heuristically tackled Tailored perturbation operators are devised and applied to both subproblems to escape from nonpromising solutions Although the topdown solution scheme is not new and has already been used by many researchers as will be shown in the literature review please see Section 2 the new major feature introduced here besides the introduction of the richer variant of the PRP is the adoption of an implicit cost to estimate delivery decisions This cost acts as a guide to obtain improved solutions by estimating how each delivery will affect the current costs of production storage and routing costs Moreover we devised different perturbation mechanisms prioritizing either tactical or operational decisions or both Finally we carried out extensive computational experiments to assess the proposed al gorithms The results show that the methods outperformed CPLEX in finding more better and faster solutions having the algorithm with the perturbation procedure that prioritizes operational level decisions outperformed the others The remainder of this paper is organized as follows Section 2 reviews recent works on PRP while Section 3 provides a formal de scription of the proposed problem and the adopted notation and defini tions The devised matheuristic approaches are presented in Section 4 whereas Section 5 presents the computational results and extensive statistical analysis Section 6 concludes with our final remarks 2 Literature review According to Bard and Nananukul 2010 PRP has four critical decisions to be made a how many items to manufacture each period b when to visit each customer c how much to deliver to a customer during a visit and d the delivery routes to use We can also include as a fifth decision e the amount of each product to be stored at the manufacturer plants and clients These decisions belong to two different decision levels a tactical responsible for production storage and backordering plans and an operational one responsible for prod uct distribution and routing The work of Adulyasak et al 2015b extensively studied PRP by listing the main devised solution methods while the works of Reimann et al 2014 and DíazMadroñero et al 2015 have discussed PRP and its many variants Here we discuss the main works that came after these literature reviews Adulyasak et al 2014b were the first to apply the adaptive large neighborhood search ALNS to the classic PRP achieving results that outperformed the existing heuristics at that time for most tested in stances BeloFilho et al 2015 addressed a case with perishable goods They also adopted ALNS but with a simulated annealing acceptance criterion Decomposition approaches are among the most intuitive and widely adopted solution strategies in the PRP literature Adulyasak et al 2015b Following this concept Absi Archetti DauzèrePérès and Feillet 2014 designed a twophase iterative mixedinteger program MIPbased heuristic to solve the uncapacitated version of the prob lem When compared with Archetti et al 2011 and Adulyasak et al 2014b their algorithm outperformed both However the Absi et al 2014 procedure was later outperformed in most sets of solution qual ity but spending more time by Solyalı and Süral 2017 while Russell 2017 achieved equivalent or superior results Working with a five phase heuristic Solyalı and Süral 2017 considered a tour idea solved via a TSP Russell 2017 adopted a matheuristic based on predeter mined or seeded routes They used the concept of a reactive tabu search to improve the previously generated routes with artificial demands Miranda et al 2018 and Qiu Wang Xu et al 2018 split the classical PRP into lotsizing and vehiclerouting problems wherein the former was solved with a MIP solver whereas the latter heuristically Miranda et al 2018 presented a formulation with products sharing the production line during some periods Their method achieved near optimal solutions when compared with CPLEX Qiu Wang Xu et al 2018 applied a skewed general variable neighborhood search and guided variable neighborhood descent to the delivery and routing decisions respectively Their method outperformed previous existing methods for solving the benchmark instances of Archetti Bertazzi Laporte and Speranza 2007 standard and high transportation costs and Boudia et al 2007 large size with comparable computational times Brahimi and Aouam 2016 were the first authors to address the PRP considering the possibility of backorders They used a relaxand fix heuristic and used it to build a solution for the lotsizing problem Their heuristic outperformed the CPLEX solving the two formulations proposed while providing good solutions in a few minutes Watanabe et al 2017 used the relaxandfix heuristic to obtain an initial solution and a fixandoptimize approach to improve solutions of a classical PRP Promising results were observed for instances with up to 30 customers NevesMoreira AlmadaLobo Cordeau Guimarães and Jans 2019 solved a large multiperishableitem PRP considering time windows for a real meat producer with a threephase fixandoptimize methodology Adulyasak et al 2014a introduced multivehicle formulations with and without vehicleindex for the classical PRP They solved the formu lations using a branchandcut BC algorithm The vehicleindexed formulation was superior in obtaining optimal solutions whereas the nonvehicleindex formulation generally provided better lower bounds on larger instances that did not reach optimality The BC algorithm attained better lower bounds at the root node of the branchandbound search tree as well as the best overall gaps Adulyasak et al 2015a were the first authors to solve the PRP under demand uncertainty considering scenarios of two or multiple stages They devised a Benders 1962 decomposition algorithm embedded in a BC scheme The computational results showed that their method was more efficient in handling instances with a greater number of scenarios than CPLEX Recently there has been concern about the impact of supply chain activities on the environment Carbon footprints and greenreverse logistics have attracted the attention of several researchers A PRP Computers Industrial Engineering 171 2022 108369 3 A Fortes et al approach for carbon and cap trade and partial delivery with lost sales was proposed by Qiu et al 2017 They assumed an objective function wherein the costs of the generated emissions were proportional to the carbon price Their model proved capable of simultaneously reducing the emission levels of carbon dioxide and operational costs Fang et al 2017 investigated a PRP with a reverse logistic also considering simultaneous pickup and delivery assumption to reduce carbon emis sions Qiu Ni et al 2018 studied a PRP with reverse logistics and remanufacturing conditions In recent years more variants of PRP have been proposed that consider aspects to make the problem more appropriate to certain situations Chitsaz Cordeau and Jans 2019 addressed the assembly routing problem that combines the production decisions with inbound logistics activities The proposed algorithm was a threephase decompo sition matheuristic that outperformed the stateoftheart heuristics at that time Considering a multiproduct PRP with outsourcing Li Chu Chu and Zhu 2019 developed a threelevel heuristic The efficiency of the algorithm was proved using tests over large instances and providing new best solutions for 283 of the 1530 benchmark instances Avci and Yildiz 2019 addressed the PRP with visit spacing policy to obtain higher quality service solutions To solve this problem a matheuristic approach based on a multistart with an MIP iterative local search was proposed When solving some standard PRP instances their algorithm performed better than those of Absi et al 2014 Adulyasak et al 2014b Li et al 2019 Qiu Wang Xu et al 2018 Solyalı and Süral 2017 and Chitsaz et al 2019 Dayarian and Desaulniers 2019 proposed a caterer PRP that con sidered shortlifespan products that are produced and consumed within a day A branchpriceandcut BPC algorithm was devel oped to solve small instances of the problem They also devised three matheuristics capable of solving larger instances nonetheless all had similar average gaps Avci and Yildiz 2020 addressed the first PRP considering trans shipment with single or multiple vehicles The problem was solved using a fast mathematical programmingbased heuristic that obtained highquality solutions Golsefidi and Jokar 2020 developed a robust approach to a PRP with uncertainty of the demand as well as simul taneous pickup and delivery The problem was solved with simulated annealing and genetic algorithms performing the latter better than the former Li Chu Côté Coelho and Chu 2020 introduced a multi plant PRP with perishable products with packaging considerations The problem was solved using a hybrid method and BC algorithms with the hybrid method obtaining good solutions within a short compu tational time Motivated by a reallife problem in the petrochemical industry Schenekemberg Scarpin Pecora Jr Guimarães and Coelho 2021 presented a twoechelon PRP An exact algorithm wherein local search procedures were incorporated within a BC framework was proposed to solve this problem Table 1 summarizes the information of reviewed papers with re gards to the characteristics on production inventory distribution and routing besides the solution approaches Production considers whether there are single or multiple plants products as well as capacity con straints The inventory defines whether there are limits to storage as well as whether it works with the maximum level or the orderupto policies Coelho et al 2013 Distribution considers whether the deliv ered loads can be realized by different vehicles in a given time period 𝑡 and if the demand can be backordered to a future period Routing determines whether the fleet is homogeneous or heterogeneous in terms of capacity and if the number of vehicles is unlimited limited single or multiple DíazMadroñero et al 2015 The solution distinguishes if the adopted approach is exact heuristic metaheuristic or hybrid blending exact and heuristic methods in some manner On observing Table 1 it is clear that there is a lack of studies that consider the ideas of backordering deliveries and routing decisions for a heterogeneous fleet features that are commonly observed in practice In addition hybrid methods that take advantage of combina tions of mathematical procedures and heuristics have been extensively explored In this manner our work proposes a richer problem and solution algorithms that contribute to fulfilling these gaps 3 Problem description and proposed formulation The richer PRP here addressed relies on the following notations and definitions Let 1 𝑇 be a set of finite and discrete planning horizon Let 1 𝑛 be the set of clients supplied by a plant with products of the set 1 𝑃 Each client 𝑖 demands 𝑑𝑡 𝑖𝑘 0 units of product 𝑘 in period 𝑡 Let 0 be the set of nodes with the plant labeled 0 and its copy 𝑛 1 Every time period 𝑡 the plant manufactures product 𝑘 while satisfying production capacity 𝐶𝑘 it incurs a setup 𝑙𝑘 and 𝑢𝑘 unit costs A product 𝑘 can be stored at the plant or clients up to the limit of 𝑈𝑘 𝑖 units 𝑖 but incurs a unitary holding cost ℎ𝑡 𝑖𝑘 for each period 𝑡 Every time the demand of client 𝑖 for product 𝑘 cannot be fulfilled in a timely fashion in a period 𝑡 the unfulfilled quantity can be backordered but at a penalty unitary cost 𝐵𝑘 𝑖 A limited heterogeneous fleet given by a set of vehicles 1 𝑉 delivers the products via routes that start at the plant node 0 and end at its copy 𝑛1 at the end of each period 𝑡 Every vehicle 𝑣 that starts a route in a time period 𝑡 incurs an activation cost 𝑒𝑣 and the routing cost associated with the arcs of the directed graph 0 where 0 𝑛 1 and 𝑖 𝑗 𝑖 𝑗 0 𝑖 𝑗 Every arc 𝑖 𝑗 has a routing cost 𝑐𝑖𝑗 0 and traversing time 𝑎𝑖𝑗 0 To unload the products from the vehicles at the clients it takes a total of 𝑠𝑖 units of service time To ease the reading Table 2 lists the parameters of the problem with parameters 𝑀 𝑡 𝑘 and 𝑀𝑣𝑡 𝑖𝑘 representing the maximum quantities for manufacturing and delivering the product 𝑘 to client 𝑖 using vehicle 𝑣 in period 𝑡 respectively A vehicleindexed formulation is here proposed This model extends the model presented in Adulyasak et al 2015b to our rich production routing problem RPRP It uses the following decision variables The productionrelated variables are variable 𝑦𝑡 𝑘 0 1 is equal to one if a setup is realized to manufacture 𝑝𝑡 𝑘 0 units of product 𝑘 in period 𝑡 and zero otherwise The variable 𝐼𝑡 𝑖𝑘 0 indicates the inventory level for product 𝑘 at node 𝑖 in period 𝑡 The variable 𝑏𝑡 𝑖𝑘 0 represents the quantity of unfilled demand for product 𝑘 at node 𝑖 in period 𝑡 The counterpart routing variables are variable 𝑔𝑣𝑡 0 1 is equal to one if vehicle 𝑣 is set to serve the clients demands of period 𝑡 otherwise it is zero The variable 𝑧𝑣𝑡 𝑖 0 1 is equal to one if client 𝑖 is visited by vehicle 𝑣 to deliver 𝑞𝑣𝑡 𝑖𝑘 0 units of product 𝑘 in period 𝑡 otherwise it is zero Finally the variable 𝑥𝑣𝑡 𝑖𝑗 0 1 is equal to one if arc 𝑖 𝑗 is used by vehicle 𝑣 in period 𝑡 otherwise it is zero The variable 𝑤𝑣𝑡 𝑖𝑗 0 represents the arrival time at node 𝑗 0 of vehicle 𝑣 after traversing arc 𝑖 𝑗 in period 𝑡 To facilitate the representation Table 3 lists all decision variables With the aforementioned notation and definitions the studied RPR can be formulated as min 𝑡 𝑘 𝑙𝑘𝑦𝑡 𝑘 𝑢𝑘𝑝𝑡 𝑘 𝑖 ℎ𝑘 𝑖 𝐼𝑡 𝑖𝑘 𝑖 𝐵𝑘 𝑖 𝑏𝑡 𝑖𝑘 𝑣 𝑒𝑣𝑔𝑣𝑡 𝑖𝑗 𝑐𝑖𝑗𝑥𝑣𝑡 𝑖𝑗 1 st 𝑝𝑡 𝑘 𝑀 𝑡 𝑘𝑦𝑡 𝑘 𝑘 𝑡 2 𝐼𝑡 0𝑘 𝐼𝑡1 0𝑘 𝑝𝑡 𝑘 𝑖 𝑣 𝑞𝑣𝑡 𝑖𝑘 𝑘 𝑡 3 𝐼𝑡 𝑖𝑘 𝐼𝑡1 𝑖𝑘 𝑣 𝑞𝑣𝑡 𝑖𝑘 𝑑𝑡 𝑖𝑘 𝑏𝑡 𝑖𝑘 𝑏𝑡1 𝑖𝑘 𝑖 𝑘 𝑡 4 𝐼𝑡 𝑖𝑘 𝑈𝑘 𝑖 𝑖 𝑘 𝑡 5 𝑞𝑣𝑡 𝑖𝑘 𝑀𝑣𝑡 𝑖𝑘𝑧𝑣𝑡 𝑖 𝑖 𝑘 𝑣 𝑡 6 𝑘 𝑖 𝑞𝑣𝑡 𝑖𝑘 𝑄𝑣𝑔𝑣𝑡 𝑣 𝑡 7 𝑣 𝑧𝑣𝑡 𝑖 1 𝑖 𝑡 8 Computers Industrial Engineering 171 2022 108369 4 A Fortes et al Table 1 Summary of the reviewed papers AUTHORS CHARACTERISTICS PRODUCTION INVENTORY DISTRIBUTION ROUTING SOLUTION Plants Product Capacity Policy Capacity Partial Backordering Vehicles Fleet Capacity Type Approach Absi et al 2014 S S no ML yes no no M HO yes HR Decomposition Adulyasak et al 2014a S S yes MLOU yes no no M HO yes EX BC Adulyasak et al 2014b S S yes ML yes no no M HO yes HR ALNS Adulyasak et al 2015a S S yes MLOU yes no no M HO yes EX BendersBC BeloFilho et al 2015 M M yes yes no no M HO yes HR ALNS Brahimi and Aouam 2016 S M yes ML yes no yes L HO yes HR RF Fang et al 2017 S S yes ML yes no no S HO yes EX BC Qiu et al 2017 S S yes ML yes yes no M HO yes EX BP Russell 2017 S S yes ML yes no no M HO yes HY MIP metaheuristic Solyalı and Süral 2017 S S yes ML yes no no M HO yes HR Decomposition Watanabe et al 2017 S M yes ML yes no no L HO yes HR RFFO Darvish Archetti and Coelho 2018 S S yes ML yes no no S HO yes EX BC Miranda et al 2018 S M yes ML yes no no M HE yes HY Decomposition Qiu Ni et al 2018 M S yes ML yes no no M HO yes EX BC Qiu Wang Fang et al 2018 S S yes ML yes no no M HO yes EX BC Qiu Wang Xu et al 2018 S S yes ML yes no no M HO yes HY VNSVND Avci and Yildiz 2019 S S yes ML yes no no SM HO yes HY Multistart and ILS Chitsaz et al 2019 S S yes ML yes no no M HO yes HY Threephase matheuristic Dayarian and Desaulniers 2019 S M yes ML yes no no M HO yes EXHY BPCMatheuristics Li et al 2019 S M yes ML yes no no M HO yes HR Threelevel heuristic NevesMoreira et al 2019 S M yes ML yes no no M HE yes HR FO Avci and Yildiz 2020 S S yes ML yes no no S HO yes HY Mathematical based heuristic Golsefidi and Jokar 2020 S S yes ML yes no no U HO yes HR GASA Li et al 2020 M S yes ML yes no no M HO yes HYEX BCMatheuristic Schenekemberg et al 2021 M M yes MLOU yes no no M HO yes EXHY BCLS This study S M yes ML yes no yes L HE yes HY ILS matheuristic S single M multiple ML maximum level OU orderupto HO homogeneous HE heterogeneous L limited U unlimited EX exact HR heuristic or metaheuristic HY hybrid Table 2 Problem parameters Notation Meaning 𝐶𝑘 𝑙𝑘 𝑢𝑘 Production capacity setup and unitary costs of product 𝑘 𝑈 𝑘 𝑖 ℎ𝑘 𝑖 𝐵𝑘 𝑖 Inventory capacity holding and backorder costs of product 𝑘 for client 𝑖 𝑒𝑣 𝑄𝑣 Cost and capacity of each vehicle 𝑣 𝑐𝑖𝑗 𝑎𝑖𝑗 Connection costs and traveling times of each arc 𝑖 𝑗 𝐻 𝑠𝑖 Maximum riding time and service times of each customer 𝑖 𝑑𝑡 𝑖𝑘 Periodic demand of each client 𝑖 of product 𝑘 in period 𝑡 𝑀 𝑡 𝑘 min𝐶𝑘 𝑖 𝑇 𝑒𝑡 𝑑𝑒 𝑖𝑘 𝑡 𝑘 𝑀𝑣𝑡 𝑖𝑘 min𝑈 𝑘 𝑖 𝑄𝑣 𝑇 𝑒𝑡 𝑑𝑒 𝑖𝑘 𝑖 𝑘 𝑣 𝑡 Table 3 Decision variables of the vehicleindexed formulation Notation Meaning 𝑏𝑡 𝑖𝑘 0 Quantity backordered at customer 𝑖 of product 𝑘 in period 𝑡 𝐼𝑡 𝑖𝑘 0 Inventory level at node 𝑖 of product 𝑘 at the end of period 𝑡 𝑝𝑡 𝑘 0 Quantity of product 𝑘 manufactured in period 𝑡 𝑞𝑣𝑡 𝑖𝑘 0 Quantity of product 𝑘 delivered to customer 𝑖 in period 𝑡 using vehicle 𝑣 𝑤𝑣𝑡 𝑖𝑗 0 The arrival time at node 𝑗 0 after vehicle 𝑣 traverses arc 𝑖 𝑗 in period 𝑡 𝑔𝑣𝑡 0 1 It is equal to one if vehicle 𝑣 performs a route in period 𝑡 otherwise it is zero 𝑥𝑣𝑡 𝑖𝑗 0 1 It is equal to one if arc 𝑖 𝑗 is used by vehicle 𝑣 in period 𝑡 otherwise it is zero 𝑦𝑡 𝑘 0 1 It is equal to one if the product 𝑘 is manufactured in period 𝑡 otherwise it is zero 𝑧𝑣𝑡 𝑖 0 1 It is equal to one if node 𝑖 is visited by vehicle 𝑣 in period 𝑡 otherwise it is zero 𝑧𝑣𝑡 0 𝑧𝑣𝑡 𝑛1 2𝑔𝑣𝑡 𝑣 𝑡 9 𝑖 𝑥𝑣𝑡 0𝑖 𝑧𝑣𝑡 0 𝑣 𝑡 10 𝑖 𝑥𝑣𝑡 𝑖𝑛1 𝑧𝑣𝑡 𝑛1 𝑣 𝑡 11 𝑖𝑗 𝑥𝑣𝑡 𝑖𝑗 𝑧𝑣𝑡 𝑖 𝑖 𝑣 𝑡 12 𝑗𝑖 𝑥𝑣𝑡 𝑗𝑖 𝑧𝑣𝑡 𝑖 𝑖 𝑣 𝑡 13 𝑤𝑣𝑡 0𝑖 𝑎0𝑖𝑥𝑣𝑡 0𝑖 𝑖 𝑣 𝑡 14 𝑖𝑗 𝑤𝑣𝑡 𝑖𝑗 𝑗𝑖 𝑤𝑣𝑡 𝑗𝑖 𝑖𝑗 𝑠𝑖 𝑎𝑖𝑗𝑥𝑣𝑡 𝑖𝑗 𝑖 𝑣 𝑡 15 0 𝑤𝑣𝑡 𝑖𝑗 𝐻𝑥𝑣𝑡 𝑖𝑗 𝑖 𝑗 𝑣 𝑡 16 Computers Industrial Engineering 171 2022 108369 5 A Fortes et al 𝑝𝑡 𝑘 𝐼𝑡 𝑖𝑘 𝑏𝑡 𝑖𝑘 𝑞𝑣𝑡 𝑖𝑘 0 𝑦𝑡 𝑘 0 1 𝑖 𝑘 𝑣 𝑡 17 𝑔𝑣𝑡 𝑧𝑣𝑡 𝑖 0 1 𝑖 𝑣 𝑡 18 𝑥𝑣𝑡 𝑖𝑗 0 1 𝑤𝑣𝑡 𝑖𝑗 0 𝑖 𝑗 𝑣 𝑡 19 The objective function 1 minimizes the production setup hold ing backordering fleet and routing costs Constraints 2 limit the quantity to be produced for each item either by the plants production capacity for that item or by its total remaining demand for horizon planning Constraints 35 control the inventory level of each prod uct in each period For each product 𝑘 both constraints 3 and 4 guarantee the balance of the inventory level at the end of period 𝑡 at the plant and customers respectively The maximum inventory levels at the factory and clients are stated via 5 The distribution and routing constraints are given by 616 The quantity of each item 𝑘 delivered to each client 𝑖 in each period 𝑡 is limited by constraints 6 Constraints 7 assure that the vehicles capacity are not exceeded Constraint 8 limits the maximum number of times that a client 𝑖 is visited per period 𝑡 whereas constraint 9 stipulates that if a vehicle 𝑣 leaves the plant node 0 it must return to its node copy 𝑛 1 Constraints 1013 are the degreerelated constraints for the plant its node copy and the clients visited by the same vehicle 𝑣 in a period 𝑡 Proposed by Bianchessi et al 2018 constraints 1416 to bind the riding time of the vehicle while eliminating the subtours The vehicle arrival times at the clients after leaving the plant are properly determined by constraints 14 while constraints 15 balance the ve hicles time continuity Constraints 16 bound the vehicles maximum riding time Constraints 1719 are related to the variables domain We assume that the backorders in the last period are equal to zero as in Brahimi and Aouam 2016 Herein formulation 119 is named as VINDX for sake of simplicity 4 Efficient matheuristics approaches As the studied RPRP inherits the computational complexity of the embedded capacity routing problem within it that is the studied RPRP is NPhard we herein present the three devised hybrid algorithms to solve it Hybrid methods seek to blend the speed of heuristics with the determinism of the exact methods Such a solution framework is normally called a matheuristic Matheuristic approaches explore math ematical programming techniques in metaheuristic frameworks or on granting to mathematical programming approaches the crossproblem robustness and constrainedCPUtime effectiveness which character ize metaheuristics Caserta Voß 2009 The devised algorithms work in twolevel decisionmaking within an iterated local search ILS framework Lourenço Martin Stützle 2003 The tactical level is re sponsible for the production and inventory plans while the operational level plans the distribution of the products and routing of the vehicles The three proposed hybrid methods follow a topdown hierarchical approach 41 Building a solution Algorithm 1 provides an initial solution based on the problem information and parameter to the topdown methods It constructs a feasible plan for the tactical problem and then proposes delivery routes for each period Throughout the execution of the algorithm the parameter 𝛥𝑣𝑡 𝑖𝑘 𝑖 𝑘 𝑣 𝑡 is a vector that has the delivery estimated costs which plays an important role in our methods 𝛥𝑣𝑡 𝑖𝑘 estimates the impact that one unit of the delivery lot 𝑞𝑣𝑡 𝑖𝑘 will have on the production inventory and transportation decisions To obtain an initial solution to the problem Algorithm 1 computes 𝛥𝑣𝑡 𝑖𝑘 as shown in Eq 20 The visitation cost 𝛥𝑖 is calculated as in Qiu Wang Xu et al 2018 The other components consider the cost of transporting a unit of product 𝑘 by each vehicle 𝑣 using the fraction 𝑒𝑣 𝑄𝑣 and if produced whether this unit will be stored either at the plant or at the customer The parameter 𝜔 1 maxWeight belong ing to a discrete and uniform distribution adds a random component to the problem to help escape nonpromising search areas 𝛥𝑣𝑡 𝑖𝑘 𝜔 𝛥𝑖 𝑒𝑣 𝑄𝑣 𝑢𝑘 minℎ𝑘 0 ℎ𝑘 𝑖 𝑖 𝑘 𝑣 𝑡 20 Algorithm 1 Build and optimization Data Problem information costs 𝚫 1 𝑞 BuildAndSolvePIv𝚫 2 VRt 3 for 𝑡 do 4 vr𝑞𝑡 BuildVRP𝑞𝑡 5 VRt VRt NRSvr𝑞𝑡 6 end 7 𝑠 PIv VRt 8 return 𝑠 Algorithm 1 solves in line 1 the toptier problem PIv with the indirect costs 𝛥𝑣𝑡 𝑖𝑘 as expressed in 2123 via CPLEX The variables and constraints of formulation 2123 have the same meaning as before please see Table 3 With the adoption of solving PIv without the visitation 𝑧 and vehicle activation 𝑔 binary variables is possible because carries both routing information min 𝑡 𝑘 𝑙𝑘𝑦𝑡 𝑘 𝑢𝑘𝑝𝑡 𝑘 𝑖 ℎ𝑘 𝑖 𝐼𝑡 𝑖𝑘 𝑖 𝐵𝑘 𝑖 𝑏𝑡 𝑖𝑘 𝑣 𝛥𝑣𝑡 𝑖𝑘𝑞𝑣𝑡 𝑖𝑘 21 st 25 𝑘 𝑖 𝑞𝑣𝑡 𝑖𝑘 𝑄𝑣 𝑣 𝑡 22 0 𝑞𝑣𝑡 𝑖𝑘 𝑀𝑣𝑡 𝑖𝑘 𝑖 𝑘 𝑣 𝑡 23 After solving PIv the attained loads are summed into parameter 𝑞𝑡 𝑖 𝑘 𝑣 𝑞𝑣𝑡 𝑖𝑘 representing the total amount of mixed product to be delivered to client 𝑖 in period 𝑡 Next the VRt problem is initialized Alg 1 line 2 This corresponds to the set of independent routing solutions for periods 𝑡 For each period the BuildVRP function generates feasible routes vr𝑞𝑡 Alg 1 line 4 This is managed via randomly selected con structive methods chosen within the following the ClarkeWright sav ings heuristic parallel or sequential see Clarke Wright 1964 and Lysgaard 1997 sequentiallexicographic and loadordered insertions For all these methods the vehicles with the largest capacities are filled first The reason for using these four constructive methods is to allow diversification in the search space In line 5 the routing solution vr𝑞𝑡 is improved by the procedure neighborhood routing search NRS procedure which is described later in Section 42 with the resulting routing solution added to the VRt list of routes for each period Although the first two constructive routing methods are well known in the literature the last two are proposed herein The sequentiallexicographic insertion starts from a node defined as 1 and attempts to add it to a vehicle 𝑣 considering the feasibility of vehicle capacity and riding time The procedure investigates all customers until no further addition to 𝑣 is possible When a vehicle is fully loaded or reaches the maximum riding time a new vehicle is activated The procedure stops after the allocation of all clients The largest vehicles are activated first The loadordered insertion algorithm takes the individual lots 𝑞𝑖 of each node 𝑖 and orders them from the largest to the smallest load The loads are then iteratively added to the largest vehicle available until no further additions can be performed Computers Industrial Engineering 171 2022 108369 6 A Fortes et al A new vehicle is activated and the process repeated until all loads are loaded For instance suppose that there are five retailers i1 i5 with delivery lots equal to q1 3 q2 2 q3 4 q4 2 q5 1 re spectively and that there are two vehicles v1 v2 with capacities equal to Q1 7 Q2 5 The sequentiallexicographic solution is the vehicle v1 visiting retailers i1 i2 and i4 while vehicle v2 visits i3 and i5 For the loadordered insertion procedure the attained solution is a vehicle v1 visiting retailers i3 and i1 while vehicle v2 visits i2 i4 and i5 42 Neighborhood routing search NRS The NRS procedure shown in Algorithm 2 applies a set of local searches to a given routing solution It is based on the random variable neighborhood descent procedure Penna Subramanian Ochi 2013 Souza Coelho Ribas Santos Merschmann 2010 that prevents method favoritism and applies each routing local search as a neighbor hood structure The adopted routing local searches are divided into two groups inter and intraroutes Kytöjoki Nuortio Bräysy Gendreau 2007 First we set RS for interroute searches line 1 While this set RS is not empty that is while solution improvements are obtained lines 210 an interroute local search rs is randomly chosen from the set RS and applied to the current solution vr𝑞𝑡 line 4 In the case of improvement an intraroute search procedure was performed line 6 Otherwise the selected local search rs is removed from the set RS line 8 Algorithm 2 Neighborhood routing search NRS Data vehiclerouting solution vr𝑞 1 RS interroute local search list 2 while RS do 3 rs 𝑟𝑎𝑛𝑑𝑜𝑚RS 4 vrq InterRouters vrq 5 if 𝑓vrq 𝑓vrq then 6 vrq IntraRoutevrq 7 else 8 RS RS rs 9 end 10 end 11 return vrq Seven neighborhood moves of interroute local searches were adopted They are Penna et al 2013 as follows Shift10 Swap11 Shift20 Swap21 Swap22 Cross and KShift K 345 These moves are realized only if a new route arrangement is feasible and improves the routing cost In the shift 10 move customer 𝑖 is transferred from route 𝑟1 to route 𝑟2 The move Swap11 picks a client 𝑖 from route 𝑟1 and client 𝑗 from route 𝑟2 and exchanges them In Shift20 move two adjacent customers 𝑖 𝑗 from route 𝑟1 are transferred to route 𝑟2 The move Swap21 exchanges two adjacent customers 𝑖 𝑗 from route 𝑟1 by customer 𝑘 from route 𝑟2 The move Swap22 exchanges a pair 𝑖 𝑗 of adjacent customers belonging to route 𝑟1 by a second pair of adjacent customers 𝑘 𝑙 from route 𝑟2 The cross move removes arcs connecting adjacent 𝑖 𝑗 clients belonging to a route 𝑟1 and the adjacent clients 𝑘𝑙 belonging to route 𝑟2 Then the arcs connecting 𝑖 𝑙 and 𝑘 𝑗 are inserted The move KShift removes a subset of consecutive 𝐾 customers from route 𝑟1 and is transferred to the end of route 𝑟2 It should be pointed out that this move is also applied if the second route is empty Intraroute local searches work reordering the visited nodes by each route We adopted five of them 1 2 or 3point move 2opt and Oropt Or 345 Groër et al 2010 described them as follows the 1point move relocates an existing node into a new position The 2point move swaps the position of two nodes The 3point move swaps the position of a pair of adjacent nodes with the position of a third node The 2opt move removes two edges from the solution and replaces them with two new edges The Oropt move removes a string of nodes and inserts it into a new position To achieve better results we also applied the very largescale neigh borhood search VLNS1 proposed by Ahuja et al 2000 It replaces both moves Shift10 and Swap11 with substantially more neigh bors This occurs because VLNS simultaneously contemplates swap and shift moves not only between a pair of routes as done by the regular shift 10 and swap 11 moves but among all sets of routes if the move is feasible Its design is extremely efficient for the family of set partitioning problems When the VLNS neighborhood structure can be classified as a networkflowbased improvement algorithm built on an improvement graph that contains all feasible moves of retailers between routes The improvement graph for a neighborhood with multiple ex changes is defined over a feasible routing solution 𝑠 for the problem which is represented by 𝐆𝑠 Let 𝑟𝑖𝑗 be a route that contains the client 𝑖𝑗 The graph 𝐆𝑠 is a directed graph with 𝑛𝑉 1 nodes where 𝑛 is the number of clients belonging to the 𝑉 routes and one additional artificial node In graph 𝐆𝑠 each node 𝑖1 𝑖𝑛 corresponds to the customers while nodes 𝑖𝑛1 𝑖𝑛𝑚 correspond to vehiclesroutes The artificial node is given by 𝑖𝑛𝑚1 A directed arc 𝑙 𝑘 𝐆𝑠 indicates that client 𝑖𝑙 leaves its current route and is transferred to the route that contains customer 𝑖𝑘 that is to route 𝑟𝑖𝑘 Simultaneously the client 𝑖𝑘 leaves 𝑟𝑖𝑘 To build 𝐆𝑠 all pairs of elements 𝑖𝑙 𝑖𝑘 𝑠 are considered An arc 𝑙 𝑘 is added to 𝐆𝑠 if and only if i the clients 𝑖𝑙 and 𝑖𝑘 belong to different routes and ii the route 𝑟𝑖𝑘𝑖𝑘𝑖𝑙 is feasible The cost 𝑐𝑙𝑘 of the arc 𝑙 𝑘 is defined by the difference in the cost of the modified route and its previous version that is 𝑐𝑟𝑖𝑘 𝑖𝑘 𝑖𝑙 𝑐𝑟𝑖𝑘 Let 𝐷 be a directed cycle on the improvement graph 𝐆𝑠 if the elements that compose 𝐷 belong to different partitions A valid cycle is defined as a directed cycle with a negative cost of 𝐆𝑠 Thus a valid cycle corresponds to a cyclic or path exchange leading to an improvement in the objective function of the problem concerning all constraints If no improvement is found the search stops To achieve these improvements efficiently identifying these valid cycles in 𝐆𝑠 is necessary Here it is performed with a labelcorrecting algorithm that obtains the minimum path from a given artificial source node to all nodes of the network It is possible that valid cycles in 𝐆𝑠 may not be found nonetheless a valid path exchange can be found Cyclic and path exchanges are a generalization of the twoexchange neighborhood Thompson Orlin 1989 Given a subset of routes the cyclic exchange removes one node from each route that composes this subset and exchanges them Fig 1 illustrates the cyclic and path exchanges Shifting the nodes 𝑖𝐴 𝑖𝐷 𝑖𝐺 and 𝑖𝐽 among their respective partitions 𝑃1 𝑃2 𝑃4 and 𝑃3 in the sequence 𝑖𝐴 𝑖𝐷 𝑖𝐽 𝑖𝐺 𝑖𝐴 characterizes a cycle exchange Fig 1a While moving the nodes 𝑖𝐴 𝑖𝐷 and 𝑖𝐺 in the sequence 𝑖𝐴 𝑖𝐷 𝑖𝐺 among partitions 𝑃1 𝑃2 and 𝑃3 characterizes a path exchange Fig 1b The cardinality of the routes involved in the exchange is maintained for the cyclic but for the path one route has its cardinality increased by one while the others have it decreased Further details regarding the VLNS method can be obtained in Ahuja 2017 Ahuja et al 1998 2000 43 Perturbation operators Local searches can become trapped in basins of attraction leading to nonpromising neighborhoods To avoid this a solution perturbation is applied and perturbation operators must be strong enough to allow the exploration of new regions and to prevent their efforts to be undone by 1 The VLNS algorithm must not be confused with the large neighborhood search LNS proposed by Shaw 1998 LNS belongs to the VLNS class of heuristics Pisinger Ropke 2010 Computers Industrial Engineering 171 2022 108369 7 A Fortes et al Fig 1 Very largescale neighborhood exchanges a local search Lourenço Martin Stützle 2019 This study considers these recommendations and outlines three algorithms having different perturbation strategies The main idea of these mechanisms is to modify either or both decision levels while verifying how these changes lead to better solutions 431 Production plan operator It is possible to find in the literature diversification moves such as the ones proposed by Armentano et al 2011 and Qiu Wang Xu et al 2018 that only change the distribution plan Nonetheless to the best of our knowledge our production plan PP perturbation move is the first one that simultaneously changes the production holding and distribution plans The PP operator shifts part of the lot produced at period 𝑡 to a period 𝑡 in which the production also takes place We consider the periods wherein the production occurs but observe the following conditions i period 𝑡 corresponds to the beginning of the adjacent production stages ii there is idle production capacity at 𝑡 and iii storage capacity available at the plant to store the excess production To understand the concept of production stages let us see the fol lowing example Suppose that there are 14 periods and that production occurs at periods 𝑡 1 5 10 and 13 This implies that the production lot manufactured in 𝑡 1 fulfills the demands of subsequent periods 1 2 3 and 4 The production of period 𝑡 5 fulfills the demands of periods 5 6 7 8 and 9 The same reasoning is applied to periods of 10 and 13 Let us consider that part of the demand 𝑑𝑡 of period 𝑡 4 is back ordered in 𝑏𝑡 units to period 𝑡 5 Then 𝑑𝑡 units will be made in 𝑡 5 Let us also consider that the manufactured lot in period 𝑡 10 contemplates part of the demand of period 𝑡 13 For the sake of the example let us assume that 𝐼12 𝑑12 Thus part of the production lot of period 𝑡 10 can be transferred to 𝑡 13 to meet demand in due time that is 𝐼12 0 𝑝13 𝑑13 𝑑14 With this example it is possible to see that two possible moves can be realized on the production lots shifting the lot partially or completely to other adjacent periods where production takes place If 𝑡 𝑡 then the move anticipates production otherwise the production is delayed Eq 24 calculates the possible shifted amount which defines the feasibility of the move For each product 𝑘 it is possible to shift from period 𝑡 to an adjacent period 𝑡 the minimum value among its entire production lot 𝑝𝑡 𝑘 either the residual production capacity 𝑀 𝑡 𝑘 𝑝𝑡 𝑘 of 𝑡 or the residual holding capacity 𝑈𝑘 0 𝐼𝑡 0𝑘 at the plant during 𝑡 𝑚𝑡𝑡 𝑘 min𝑝𝑡 𝑘 𝑀 𝑡 𝑘 𝑝𝑡 𝑘 𝑈𝑘 0 𝐼𝑡 0𝑘 𝑘 𝑡 𝑡 𝑡 𝑡 24 Every time the operator PP is selected it lists all feasible moves of the production lots Subsequently it randomly draws up to 𝑚𝑎𝑥𝑃𝑒𝑟𝑡𝑃𝑃 moves for different products and periods To maintain solution feasibil ity it allows the modification of the production plan of a product 𝑘 up to one time per perturbation call The selected moves are applied to PIv by changing the bounds of the variables while considering the value of 𝛥𝑣𝑡 𝑖𝑘 432 Routing operators Exchanges on the routing plans characterize the perturbation on the operational level Five operators are responsible to modify the routes For each of them between 25 and 75 of the periods are randomly selected and perturbed but the solutions are always feasible They are separated into intraroute and interroute operators After the execution of routing perturbations the NRS Section 42 procedure was applied to the perturbed solutions The intraroutes perturbations modify one route per time We de vised the reverseroute RR lexicographic LX and randomized route RD Operator RR inverts the sense of selected routes giving more opportunity to the last visited nodes to be explored first by routing local searches It reverts the sense of up to 𝑚𝑎𝑥𝑉 routes for a period of time 𝑡 This parameter is randomly selected from the interval 0 𝑉 every time the operator is called As an example suppose that route 𝑟 visits a subset of customers in the following order 3 4 5 1 2 According to the selected intraroute operator 𝑟 would turn into 2 1 5 4 3 ap plying RR and operator LX would reorder 𝑟 in the sequence 1 2 3 4 5 while RD could generate a new 𝑟 in the order 3 1 4 2 5 Both operators LX and RD were calibrated They are executed in each selected period 𝑡 up to 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝐿𝑋 and 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝐿𝑋 times respectively The interroute perturbation operators modify two routes The split route SR operator looks for an empty vehicle and allocates to some or all retailers visited by another route For example a route that visits the following customers 1 2 4 5 3 could be divided into 1 2 4 and 53 This operator searches for an empty vehicle 𝑣1 and randomly selects another vehicle 𝑣2 𝑣1 𝑣2 which is active It then transfers as much as possible to the load from 𝑣𝑡 to 𝑣1 The KK operator randomly selects two routes and strings of sequen tial nodes with sizes 𝐾1 and 𝐾2 belonging to them These strings are exchanged between routes allocating them at the end of each route For example suppose that a route visiting clients 1 2 3 4 and a second route visiting 7 6 8 5 By applying operator KK they could be turn into 1 2 8 5 and 7 6 3 4 This perturbation operator is applied up to 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝐾𝐾 over each selected time period 𝑡 Computers Industrial Engineering 171 2022 108369 8 A Fortes et al 44 Updating Every time a perturbation operator is applied to the solution the indirect costs are updated according to Eq 25 They offer much more accurate information about the problem and assess the impact that the transported lots have on the production storage and routing components 𝛥𝑣𝑡 𝑖𝑘 𝜔 𝛥𝑖 𝛥𝑟 𝑒𝑣 𝑄𝑣 ℎ𝑘 𝑖 𝐼 𝑡 𝑖𝑘 𝐼 𝑡1 𝑖𝑘 𝑢𝑘 𝑝𝑡 𝑘 𝑝𝑡1 𝑘 𝑖 𝑘 𝑣 𝑡 25 Eq 25 is the sum of the following components The parameter 𝛥𝑖 𝑐ℎ𝑖 𝑐𝑖𝑗 𝑖 𝑗 𝑘 𝑟 that is it is equal to the sum of the costs of the corresponding arcs that connect node 𝑖 within route 𝑟 whose total length is equal to 𝛥𝑟 The second component is the fraction 𝑒𝑣 𝑄𝑣 This shows the cost of transporting a unit by each vehicle 𝑣 The component that divides the holding cost ℎ𝑘 𝑖 by the difference between the inventory levels in periods 𝑡 and 𝑡 1 estimates the impact of each load unity on the inventory from one period to the next The reasoning is analogous for the last component considering production costs and manufactured lots By adding more information to the decisionmaking at the tactical level is enriched and is able to achieve better solutions when compared with the one proposed by Qiu Wang Xu et al 2018 see Section 54 of the computational experiments 45 Iterated local search framework Our solution strategy for the studied RPRP is based on an ILS meta heuristic whose core concept involves improving a current solution by generating new initial solutions through perturbations combined with local searches Lourenço et al 2019 The devised hybrid approaches follow a topdown decisionmaking process usually adopted by production planning systems Bitran Tirupati 1993 Algorithm 3 shows the main steps of the devised ILS framework An initial solution 𝑠 is provided by Algorithm 1 and is declared as the current best solution line 1 Each iteration of the outer loop lines 220 attempts to improve 𝑠 The setup variables 𝑦 have fixed values for each period to the corresponding products being manufactured in solution 𝑠 line 4 It aims to make the resulting problem easier and quicker to be solved Line 5 initializes the set P of perturbation operators 𝜖𝜅 In the inner loop lines 619 a perturbation operator is randomly selected to avoid favoritism and applied to the current solution The chosen perturbation operator modifies either the tactical or operational part of the solution or both Recall the operators were presented in Section 43 The parameter and its updating function are at the core of the proposed ILS algorithms As previously mentioned parameter is a vector of indirect costs which is used for the tobedelivered loads and to guide the solution of the tactical problem PIv As the RPRP is solved in two stages parameter estimates the impact of the delivery lots on the production and routing plans It works as a feedback parameter to guide the solution of the toptier problem Function update Δ renews parameter at the problem PIv see Section 44 In line 10 the problem PIv is optimized considering the renewed and examined if it found a better objective function value If a better solution is obtained the search continues from it Otherwise the objective function value of the current solution is tested and if it is less than 1 𝛼 times the objective function value this solution is explored in the next iteration This step allows for diversification exploring nonlocal optimal solutions and is inspired by the skewed variable neighborhood search Hansen Mladenović 2001 A perturbation operator is discarded from the list if a solution improvement is not attained in the current iteration This strategy of the perturbation operator can be understood as a combination of the cyclic and pipe neighborhood exchange steps Hansen Mladenović Todosijević Hanafi 2017 If there is an improvement the algorithm does not immediately return to the first perturbation operator that is it can explore a different perturbation operator from the current operator unless it is drawn again However if there is no improvement the remaining operators have the same opportunity to be selected The results show that no operator was at a disadvantage to the others If the setup variables 𝑦 at the tactical problem are fixed then they are unfixed Otherwise they were fixed It works closely related to a relaxandfix procedure Pochet Wolsey 2006 Recall that with 𝑦 fixed it is expected that the problem will be solved easily and quickly Algorithm 3 Iterated local search ILS Data Problem information costs 𝚫 1 𝑠 BuildAndOptimize 𝚫 2 for 𝑖 𝑚𝑎𝑥𝐼𝑡𝑒𝑟 do 3 𝑠 𝑠 4 fix the setup variables 𝑦 5 P 𝜖𝜅 6 while P do 7 𝜖𝜅 RandomP 8 𝑠 Perturb𝑠 𝜖𝜅 9 𝑠 UpdateΔ𝑠 10 𝑠 Opt𝑠 11 if 𝑓𝑠 𝑓𝑠 then 12 𝑠 𝑠 13 if 𝑓𝑠 𝑓𝑠 then 𝑠 𝑠 14 else 15 P P 𝜖𝜅 16 if 𝑓𝑠 1 𝛼𝑓𝑠 then 𝑠 𝑠 17 unfixfix 𝑦 18 end 19 end 20 end 21 return 𝑠 The proposed algorithms were named TILS OILS and IILS ac cording to the adopted perturbation strategy tactical operational and integrated respectively TILS adopted the production plan operator PP from Section 431 OILS selects within the five routing operators RR LX RD SR and KK as presented in Section 432 Finally IILS uses all proposed perturbation operators 5 Computational results and analysis This section reports the computational experiments performed with all analyses performed for our algorithms Here we present a summary of our results with all data available in the online resource provided by the authors 51 Instance generation To the best of our knowledge there is no set of benchmark in stances available that gather all the necessary information for our formulation and algorithms Therefore we have generated a set of instances that were created using an approach similar to the works of Coelho and Laporte 2013 and Brahimi and Aouam 2016 but our instances are not the same or equivalent to theirs nonetheless since we have addressed a much richer problem here We have developed 108 instances by varying four main aspects 1 number of retailers 𝑛 10𝑐 𝑐 2 5 2 number of periods 𝑇 5 10 15 3 number of products 𝑃 3 5 7 and 4 number of available vehicles 𝑉 7 9 11 The parameters and costs were drawn from discrete uniform dis tributions within the intervals as shown in Table 4 Two different intervals were used to characterize low and high levels of demands We Computers Industrial Engineering 171 2022 108369 9 A Fortes et al Table 4 Parameters and costs values Values Meaning 𝑥𝑖 𝑦𝑖 0 1000 𝑖 The node coordinates 𝑑𝑡 𝑖𝑘 0 25 or 𝑑𝑡 𝑖𝑘 30 55 𝑖 𝑘 𝑡 𝑃 𝑇 Clients periodic demand 𝐼0 0𝑘 100 150 𝑘 Initial inventory of the plant 𝐶𝑘 𝑛𝐴𝑘 𝐴𝑘 50 140 𝑘 Production capacity 𝑢𝑘 2 8 𝑙𝑘 104𝑢𝑘 𝑘 Production setup and unitary costs 𝑈 𝑘 𝑖 𝜅𝐹𝑘 𝐹𝑘 140 190 𝜅 2 3 4 𝑖 𝑘 Holding capacity ℎ𝑘 0 1 5 ℎ𝑘 𝑖 6 10 𝑖 𝑘 𝑃 Holding cost 𝐵𝑘 𝑖 𝛾ℎ𝑘 𝑖 𝛾 8 12 𝑖 𝑘 𝑃 Backorder costs 𝐷 𝑡 𝑘 𝑖 𝑑𝑡 𝑖𝑘 Average demand per period 𝑒𝑣 500 1000 𝑄𝑣 2𝐷 𝑇 𝑉 08 1 𝑣 Vehicle cost and capacity 𝑐𝑖𝑗 𝑟𝑜𝑢𝑛𝑑 𝑥𝑖 𝑥𝑗2 𝑦𝑖 𝑦𝑗2 𝑖 𝑗 Traveling costs 𝑎𝑖𝑗 𝑐𝑖𝑗 𝑖 𝑗 The traveling times 𝐻 6000 The maximum riding time 𝑠𝑖 50 𝑖 client service time assumed as in Brahimi and Aouam 2016 that the initial backorder was equal to zero and the inventory at the clients was equal to the clients firstperiod demand At the link httpstinyurlcomrprpinstances one can find all of the proposed instances as well as the source codes for reading them and generating new ones It also has all of the source codes for the devised algorithms and the tables reporting their attained bounds 52 Parameters adjustment and implementation details All experiments were performed on an Intel Xeon CPU E5 2687 W v3 310 GHz computer with 160 GB of RAM and running Ubuntu Linux 1804 The matheuristics were coded in C and the Concert Technology API CPLEX was used to solve production planning problems The PIv had at maximum 600 s s or 5 estimated gap to run We solved each instance with the matheuristic algorithms ten times as planned by a power of sample test with a power of 95 to detect the differences in the averages We set the stopping criterion of 7200 s for all algorithms The TILS variant of Algorithm 3 only works with the PP oper ator The inner loop lines 619 was modified to repeat at most 𝑚𝑎𝑥𝐼𝑡𝑒𝑟𝑇 𝐼𝐿𝑆 number of times Every time an improvement is obtained in line 13 the counter 𝑛𝐼𝑡𝑒𝑟 0 𝑚𝑎𝑥𝐼𝑡𝑒𝑟𝑇 𝐼𝐿𝑆 is reset to 0 otherwise it is incremented by one The parameters 𝛼 maxWeight maxIter maxPertPP maxPertRD maxPertLX and maxPertKK were all calibrated using the irace package of LópezIbáñez DuboisLacoste Cáceres Birattari and Stützle 2016 an automatic configuration tool The parameter 𝛼 of Algorithm 3 allows the method to explore more configurations without straying too far from the best current solution as done in the skewed VNS Hansen Mladenović 2001 The parameter maxWeight adds an aleatory component to the problem as it modifies the vector of costs The parameters maxIter define the maximum number of iterations of the loop between lines 2 and 20 in Algorithm 3 Parameter maxPertPP bounds the maximum number of moves that perturbation operator PP realizes The parameters maxPertRD maxPertLX maxPertKK limit the number of times that the routing perturbation operators RD LX and KK are applied over the routing solution Each algorithm was executed thousand times and the irace picked the best values for each parameter Table 5 53 Impacts of backordering and a heterogeneous fleet As shown in the literature review please see Table 1 backorder and heterogeneous fleets are not features commonly adopted in PRP studies Given that a regular PRP already packs a large number of complicating assumptions the presence of additional features hinders not only the solution of a test instance but also the development of Table 5 Adopted parameter values by the devised methods Parameter TILS OILS IILS 𝛼 10 20 15 𝑚𝑎𝑥𝑊 𝑒𝑖𝑔ℎ𝑡 5 5 5 𝑚𝑎𝑥𝐼𝑡𝑒𝑟 300 300 300 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝑃 𝑃 15 5 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝑅𝐷 7 7 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝐿𝑋 3 3 𝑚𝑎𝑥𝑃 𝑒𝑟𝑡𝐾𝐾 15 10 algorithms Nevertheless their presence bridges the applicability gap of such solution frameworks In industrial applications backorders represent a risk of losing sales or even some patrons but when judiciously allowed they pose to reduce delivery or inventory costs According to Brahimi and Aouam 2016 even in a vendormanagement inventory policy backorders can occur Some causes can be listed as limited production or stor age capacity at the plant facility prioritization of other products occurrence of delay in replenishment owing to the fleet capacity or limited storage capacity at the customers However determining the quantities that can be backordered without negatively impacting costs or substantially affecting the current production plan is a nontrivial task Moreover delivery fleets can be owned or outsourced and may not always be homogeneous differing in capacity variable and fixed costs speeds and the customers that they can access Toth Vigo 2014 To illustrate how these aspects influence the obtained solutions and their costs using the CPLEX we solved the proposed set of in stances see Section 51 considering three alternative cases a no backordering and heterogeneous fleet b neither backordering nor heterogeneous fleet and c backordering and a homogeneous fleet These alternative cases are compared with the proposed RPRP case d considering simultaneous backordering and a heterogeneous fleet Some metrics are presented in Table 6 for further discussion Table 6 presents the percentage difference among the aforemen tioned cases ab and c when compared with the proposed RPRP represented by case d A negative value means that for alternative cases a worse value was found For each of analyzed cases a b c and d when solving the proposed set of instances the CPLEX was ca pable of attaining feasible solutions for 74 70 84 and 84 respectively First it suggests that allowing for backordering the demands cases c and d eases the resolution of the problems and leads to better costs and use of resources For all tested problems and cases the same production costs were obtained once the production even if backordered must satisfy cus tomers demands When compared with the proposed RPRP all al ternative cases had the worst total cost The setup costs show that Computers Industrial Engineering 171 2022 108369 10 A Fortes et al Table 6 Solutions comparison among the cases a b and c with d Cases comparison Metrics a b c Total cost 2122 688 706 Setup cost 5218 2007 1267 Production costs 000 000 000 Holding costs 181 2570 4506 Backordering costs 000 000 294289 Fleet cost 675 085 040 Traveling cost 1668 856 232 Setups 5045 1889 1289 Activated vehicles 549 000 311 Visited customers 456 342 036 Fig 2 Comparison of the new 𝛥𝑣𝑡 𝑖𝑘 with the 𝛥𝑖 of Qiu Wang Xu et al 2018 the alternative scenarios performed much worse than the RPRP case especially without the backorder policy a is approximately 52 larger than the RPRP In contrast this was the sole case that attained slightly better holding costs Observing the backordering cost metric case c was almost 3000 worse than our proposed method As we can see the absence of both backordering and heterogeneous fleet incurs higher fleet and traveling costs as well as the total cost for all three alternative cases Furthermore the disallowance of backorder policy and the adop tion of a homogeneous fleet can lead to longer travel distances a larger number of customers visited per period an increase in the number of vehicles running and times that the setups occur These settings raise environmental concerns as highlighted by Fang et al 2017 Qiu et al 2017 and Qiu Ni et al 2018 54 Importance of Recall that our parameter 𝛥𝑣𝑡 𝑖𝑘 estimates the impact of routing and distribution decisions when creating a new production plan It helps to guide CPLEX to obtain more promising production and inventory plans Here we compare our parameter with the one introduced in Qiu Wang Xu et al 2018 To the best of our knowledge this is the first time indirect costs related to delivery variables 𝑞 are used Using both theirs and our parameter all the three outlined ILS algorithms were tested on 12 instances PX0n5p15tYv where X 2 3 4 5 and Y Fig 2 brings the attained benchmark pro file Dolan Moré 2002 where the plots only use the overall best solutions of the algorithms with respect to the adopted parameter Our adaptive indirect cost 𝛥𝑣𝑡 𝑖𝑘 solid and orange lines reaches 90 of the best solution values demonstrating its superiority over 𝛥𝑖 dashed and blue lines To draw further insights into the contribution of we performed statistical tests Table 7 Here owing to the smaller sample size and the Table 7 Statistical tests for the comparison of the new 𝛥𝑣𝑡 𝑖𝑘 with 𝛥𝑖 of Qiu Wang Xu et al 2018 a KolmogorovSmirnov tests b Wilcoxon tests 𝛥 Statistic pvalue Alternative Statistic pvalue 𝛥𝑣𝑡 𝑖𝑘 10 00 default 360 31 106 𝛥𝑖 10 00 greater 6300 15 106 absence of normality please see the results of the KolmogorovSmirnov tests of Table 7a which shows that the 𝑝value is less than 005 the nonparametric Wilcoxon test WX was applied The alternative hypothesis of the WX signedrank test verifies whether the results are statistically different that is it verifies if the 𝛥𝑣𝑡 𝑖𝑘 results are better than the 𝛥𝑖 results Table 7b shows the pvalues of the Wilcoxon tests and that the results are indeed significantly different besides showing a clear superior performance of 𝛥𝑣𝑡 𝑖𝑘 over 𝛥𝑖 because both pvalues tend to zero 55 Comparing the algorithms We have assessed the performance of the algorithms with respect to the solutions and the time spent to reach them First the best solution values attained by the algorithms were compared to the upper bounds of the VINDX formulation Hereafter the three devised algorithms TILS OILS and IILS are compared We have generated instances by varying the values of our four attributes the number of customers n the number of products p the number of periods t and the number of vehicles v We have labeled the instances according to these attributes Usually for problems with fewer attributes eg the numbers of clients and candidate facilities of an uncapacitated facility location problem it is possible to correlate the difficulty of solving an instance with the values of these attributes ie with the instance size Larger values usually imply a more difficult instance to solve Nonetheless as in the present problem when one has more attributes this correlation may not be direct As we will show next depending on how the instance data affect each other smaller instances can be as difficult to solve as larger ones We have provided a posteriori analysis of the instances to propose a typology for them The idea is to classify the instances according to the attained optimality gaps as a metric We have used the achieved lower bounds after running CPLEX for six hours and the bestknown solutions BKS found by the proposed methods to compute the optimality gaps We have then divided our instances into four ad hoc main groups 1st Instances that had optimality gaps less than 20 2nd instances that had optimality gaps between 20 and 30 3rd instances that had optimality gaps between 30 and 40 and 4th instances that had optimality gaps larger than 40 It is worth highlighting that CPLEX has been able to provide an upper bound to 75 of the 108 proposed instances after six hours of computation when solving the VINDX model ie CPLEX has been unable to find a feasible integer solution for 33 instances or close to 31 of them On the other hand all the outlined matheuristics have reached better upper bounds than the CPLEX ones for all instances with an average of fewer than 30 min for the worst case Fig 3 shows the sorted attained optimality gaps colored by groups Note that a white space delimits these groups and the green horizontal line highlights the average optimality gap of 2935 Considering the 33 instances that CPLEX has not found an integer feasible solution 23 of them belong to the third and fourth groups showing that these groups with larger optimality gaps may be seen as the ones having the more difficulttosolve instances To draw some insights Table 8 shows a statistical summary of the characteristics of each group The first secondthirdfourth group is composed of 12 44439 instances having on average 3167 320536984444 customers 425 402400356 products 625 Computers Industrial Engineering 171 2022 108369 11 A Fortes et al Fig 3 Attained optimality gaps for the proposed instances Table 8 Average attribute value by instance group Group Total n p t v OPT 1st 12 3167 425 625 917 1852 2nd 44 3205 402 750 909 2491 3rd 43 3698 400 1279 891 3420 4th 9 4444 356 1389 878 4237 7512791389 periods 917 909891878 vehicles and 1852 249134204237 of optimality gap Note that the larger the number of clients and periods the larger the optimality gaps as one might expect Handling a large number of clients and periods hinders the solution process Nonetheless on the other hand in a counterintuitive manner the fewer the number of products and available vehicles the larger the optimality gaps When more vehicles are available to handle the deliveries the easier it is to find a feasible solution One can easily accommodate the delivery orders on the fleet However when fewer vehicles are available one must finetune the assessment of the trade off among the routing cost the inventory cost and the backorder cost It is more difficult to combine mixed loads while adequately serving customers When dealing with a smaller number of products the results show that these instances generally lead to greater gaps However this is not an indicative that the problem will be easily solved when a larger number of products is used One must recall that the problem itself involves other factors such as the arrangement of the demands the capacity of the inventory level the size and capacity of the fleet and the size of the planning horizon What the results show is that these factors combined are extremely difficult to deal with explaining somehow why there are just a few works on the literature addressing them A different way to visualize the aforementioned insights is by means of a bubble chart as shown in Fig 4 This figure summarizes the impacts of each of the four attributes in the attained optimality gaps It is clear that the larger the number of periods red bubbles or the number of customers n 40 or 50 the larger the optimality gaps vertical axis However in a counterintuitive manner smaller size fleets represented by the size of the bubbles may lead to larger optimality gaps Instances with fewer products p 3 have usually larger optimality gaps The summary of the upper lower bounds best known solutions and optimality gaps can be found at the Appendix A To show the difference between the average and minimum objective function values found by the outlined methods and CPLEX we have used a benchmark profile plot Performance profiles are graphic tools that aid in evaluating and comparing the performance of a set of algo rithms to solve a set of problems over a given performance measure They have been proposed by Dolan and Moré 2002 wherein more details can be obtained The benchmark profile presented in Fig 5a is adopted as a mea sure of the number of times an algorithm reaches the minimum value found for each instance also known as the bestknown solution BKS The OILS algorithm reaches 5370 of the BKS followed by the IILS 3241 TILS 1296 and VINDX 093 In general the VINDX results were between an average distance from the BKS of approxi mately 47 with a maximum of 181 registered To conclude this comparative analysis of the objective function val ues Fig 5b illustrates the comparison of the averages obtained by the ILS methods in comparison with the best values obtained by the VINDX OILS has the best average 555 followed by IILS 343 TILS 166 and distance by VINDX 093 One can perceive the pre ponderance of the heuristics in comparison to the complete model used by the solver in finding good solutions in viable computational time for these problems once the proposed methods easily found solutions for all of them within 1200 s The behavior of the matheuristics was studied by conducting statis tical tests on the objective values and running times The parametric methodology ANOVA was chosen which even with the premise of normality of the data works very well even with asymmetric distri butions as long as the sample size is significantly large which is the case All following statistical tests were performed using the MINITAB 19 The objective function values of the algorithms are studied The null hypothesis tests if the average objective function value is equal while the alternative hypothesis tests if they are not with 𝛼 005 Equality of variances of the algorithms was assumed for the analysis As 𝛼 p𝑣𝑎𝑙𝑢𝑒 005 0295 the null hypothesis is not rejected Table 9 Thus there is sufficient statistical evidence to infer that the three algorithms returned the same average value for the objective function Computers Industrial Engineering 171 2022 108369 12 A Fortes et al Fig 4 Comparison of the attained optimality gaps with the problem attributes Fig 5 Benchmark profiles for the devised methods Table 9 ANOVA tests for objective function Source F 𝑝value of 122 0295 Table 10 Averages and standard deviation for objective function Algorithm N AVG SD TILS 1080 1031022 576228 OILS 1080 995410 548640 IILS 1080 1002700 551385 AVG average SD standard deviation The Tukey simultaneous tests with 95 confidence for the average objective function was performed Table 10 lists the averages and standard deviations for the objective function value while Table 11 presents their relative differences Once the confidence interval CI contains the value 0 the difference between the average objective values can be 0 and the average of the objective function value of these methods is not considerably different It is here important to highlight the behavior of the devised methods regarding the attained objective function values One would expect that the IILS having all the perturbation mechanisms would lead to better solutions than those found by the TILS and OILS but this was not what happened Please bear in mind that when the PP operator changes a production plan it directly affects the holding and the backorder and the distribution decisions which for their turn severely impact the routing decisions It was then observed that these cascading effects of ten led to infeasible solutions that needed fixing Unfortunately fixing these solutions came at a substantial computational cost with a further sideeffect the repaired solutions often led to nonpromising basins of attraction Given the aforementioned behavior one might suggest that the PP operator presented an almost total random behavior On the other hand the operational operators responsible for perturbing the routing decisions only produced feasible solutions ie solutions Table 11 Tukey simultaneous tests for differences of objective function averages Level Diff AVG Diff SE Diff CI of 95 Tvalue 𝑝value TILS OILS 35612 24082 20753 91977 148 0301 IILS OILS 7291 24082 49074 63656 030 0951 IILS TILS 28322 24082 84687 28043 118 0468 AVG average SE standard error Diff difference CI confidence interval Table 12 Welch tests results Source F pvalue Opt 1468 000 PIv 7472 000 VRt 7080 000 Total 13304 000 that simultaneously respected all of the problems constraints Thus concerning the attained objective function values the OILS led to better solutions than the IILS and TILS while the IILS performed in between the OILS and TILS strategies The running time was analyzed by considering four different mea surements in seconds for each algorithm They are i the time until the best solution Opt ii the cumulative time spent to solve the problems PIv and iii VRt and iv total running time The time measure ii includes the time spent by the PP while iii accounts for the time spent by the routing operators RR SR KK LX and RD Oneway ANOVA was used to detect any differences between each of these times The null hypothesis is that if the average time spent is equal and the alternative hypothesis is not with 𝛼 005 No equality of variances is assumed for the analysis The Welch test is significant at the 5 level of significance that is there is sufficient statistical evidence to infer that the three algorithms have different Opt PIv VRt and total average times Table 12 once 𝛼 𝑝value and then the null hypothesis is rejected On observing Table 13 with the respective standard deviations and averages we note that the OILS algorithm has the worst average time Computers Industrial Engineering 171 2022 108369 13 A Fortes et al Table 13 Standard deviation and average of the times Time TILS OILS IILS SD AVG SD AVG SD AVG Opt 65450 72960 46250 62170 40498 72960 PIv 119870 148230 54740 108510 36918 101850 VRt 16218 12681 25816 20487 25642 22340 Total 108220 186830 45600 132290 21915 132199 SD standard deviation AVG average Table 14 Tukey simultaneous tests for differences of time averages Time Level Diff AVG Diff SE Diff CI of 95 Tvalue 𝑝value Opt TILS OILS 1079 244 508 1650 442 0000 IILS OILS 1079 244 508 1650 442 0000 IILS TILS 00 282 659 659 000 1000 PIv TILS OILS 3973 401 3034 4911 991 0000 IILS OILS 666 201 1136 196 331 0003 IILS TILS 4639 382 5532 3745 1215 0000 VRt TILS OILS 7806 928 9978 5635 841 0000 IILS OILS 185 111 74 444 167 0215 IILS TILS 9659 923 7498 11820 1046 0000 Total TILS OILS 5454 357 4617 6290 1526 0000 IILS OILS 09 154 370 351 006 0998 IILS TILS 5463 336 6249 4677 1626 0000 only for the VRt when compared with TILS algorithm when compared with the IILS algorithm they have close total and VRt times Finally when comparing the TILS and IILS algorithms the average to the best answer Opt is identical and except for the time spent in VRt all other times in TILS are worse than IILS GamesHowell paired comparison tests were performed to infer the average time difference between the algorithms The GamesHowell simultaneous test is significant at the 5 level of significance that is there is sufficient statistical evidence to infer that the difference in the average times of the three algorithms Column AVG Diff Table 14 is different from 0 Subsequently the null hypothe sis regarding the equality of the average times of the three algorithms is rejected 𝛼 𝑝value in the following pairwise comparisons For Opt OILS is different from both TILS and IILS for PIv all algorithms are different among themselves whereas for VRt and total average times TILS is different from both OILS and IILS Finally the performance of the algorithms was examined with a multiple timetotarget MTTT plot an extension of the TTT plots Aiex Resende Ribeiro 2002 2007 and capable of setting multiple instances simultaneously Reyes Ribeiro 2018 For this the worst objective function values from the five selected instances are defined as their respective targets These instances were randomly chosen from the group of problems wherein the CPLEX did not obtain any upper bound Each method was run for 100 executions until the target was found or the limit time was reached Fig 6 presents the MTTT graphic wherein it shows that the OILS algorithm in 80 of attempts reaches the target within 400 s overlapping both the IILS and the TILS algorithms with up to 435 and 580 s respectively Once again the OILS algorithm had the best performance followed by the IILS and TILS algorithms in that order The analysis is extended to the efficiency of perturbation operators and routing local searches We consider as efficiency the number of times that the adoption of the operator or local search led to solution improvements From Table 15 it can be deduced that in the TILS variation the operator PP is not efficient leading to only 143 improvement In contrast when it works together with the routing operators in IILS its efficiency increases almost 75 times A similar behavior is observed for the routing operators in OILS and IILS their performance increases by approximately 1 when coupled with PP Fig 6 Multiple TTT plot Table 15 Percentage efficiency of the perturbation operators PP RR LX RD SR KK TILS 143 OILS 493 493 495 484 484 IILS 1072 575 584 58 583 578 Considering the routing local searches the VLNS showed itself to be extremely powerful being responsible for reaching approximately 391 of all improvements followed by the SWAP21 SWAP22 SHIFT20 KSHIFT3 KSHIFT4 CROSS and KSHIFT5 with 242 151 123 37 21 20 and 14 of the time respectively This performance obtained by the VLNS is due to its well defined networkflow structure which enables many more exchanges than the other searches The set of intraroute searches had improved movements with 2opt 258 1pointmove 207 2pointmove 203 3pointmove 170 Oropt Or 3 74 Or 4 53Or 5 35 respectively After these analyses it is clear that operational decisions can strongly affect the performance of integrated problems This corrob orates the thesis that by spending more time in the prospection of the routing component most solutions converge more quickly to better results than those found focused first on the productioninventory problem This was shown in the performance of both the OILS and IILS algorithms that used operationallevel operators This leads to a reflection on whether decisionmaking can pay more attention to operational decisions and how they influence tactical decisions 551 Benchmark instances We developed tests over a set of instances proposed by Boudia et al 2007 and Archetti et al 2011 The obtained results were compared with those of Russell 2017 Solyalı and Süral 2017 and Schenekem berg et al 2021 Each algorithm was executed five times for each instance and a brief analysis is presented as follows The first sets B1 B2 and B3 from Boudia et al 2007 totaled 90 instances The algorithms had averages above 17 of the BKS for all sets demonstrating a very low adherence to these instances The sets A2 A50 and A3 A100 from Archetti et al 2011 totaled 192 instances with multiple vehicles For the A2 A3 set the results presented a better performance with an average 479 224 worse than the BKS with the best results remaining 04 03 from the BKS Future studies will involve a better calibration of the algorithm parameters to obtain more refined solutions Computers Industrial Engineering 171 2022 108369 14 A Fortes et al 6 Conclusions Owing to the highlevel computational challenge offered by the proposed rich productionrouting problem we herein developed ap proaches that split the problem into a topdown decision manner A toptier problem solves the production inventory and distribution plans with CPLEX providing an estimate of the loads to be distributed The bottom tier heuristically routed the operational level for each period Three algorithms embedded in an ILS framework were devised The first provoked perturbations at the tactical level modifying the PP and directly affecting the inventory backorder and distribution variables The second explored a set of mechanisms to modify the routing plan changing the operational level The last algorithm works by integrating perturbations at both levels An important and common point in the implementation of the algorithms was the development of an intuitive method for calculating the costs of the delivered loads This is because if there is a need to deliver a lot this considerably affects the costs related to production inventories and transportation In this manner during each iteration of the method these indirect costs are updated with the latest information from the solution Combined with the perturbation mechanisms it was able to lead to better solutions than those obtained by the commercial solver Through computational experiments the methods were compared with the proposed formulation and among themselves The CPLEX had difficulty obtaining a solution to most problems within the available time On the contrary for those instances wherein a solution was obtained CPLEX was surpassed both in running time and objective function values by the three proposed methods The results were statistically analyzed showing the importance of our indirect cost and perturbation mechanisms especially those that modified route planning The operational perturbation operator was capable of reaching most of the BKSs with less computational effort for the tested instances For future research it may be interesting though more challeng ing to add more features to PRP Features such as allowing partial deliveries or lost sales perishable products adoption of periodic time windows with demands being served within a given range of periods heterogeneous and multiple fleets thirdparty logistic operators and proportional service time may bring theory even closer to practice which will most likely lead to even greater cost savings CRediT authorship contribution statement Allexandre Fortes Conceptualization Methodology Software For mal analysis Writing original draft Ricardo Camargo Conceptual ization Supervision Writing review editing Leandro Reis Muniz Conceptualization Formal analysis Fátima Machado de Souza Lima Conceptualization Reviewing and editing Fernanda dos Reis Cota Conceptualization Formal analysis Declaration of competing interest The authors declare that they have no known competing finan cial interests or personal relationships that could have appeared to influence the work reported in this paper Acknowledgments The authors thank CNPq for financially supporting this research The second author is funded by the grant 30885220192 Appendix A Supplementary data Supplementary material related to this article can be found online at httpsdoiorg101016jcie2022108369 References Absi N Archetti C DauzèrePérès S Feillet D 2014 A twophase iterative heuristic approach for the production routing problem Transportation Science 494 784795 Absi N Archetti C DauzèrePérès S Feillet D Speranza M G 2018 Com paring sequential and integrated approaches for the production routing problem European Journal of Operational Research ISSN 03772217 2692 633646 Adulyasak Y Cordeau JF Jans R et al 2014a Formulations and branchandcut algorithms for multivehicle production and inventory routing problems INFORMS Journal on Computing 261 103120 Adulyasak Y Cordeau JF Jans R et al 2014b Optimizationbased adaptive large neighborhood search for the production routing problem Transportation Science 481 2045 Adulyasak Y Cordeau JF Jans R et al 2015a Benders decomposition for production routing under demand uncertainty Operations Research 634 851867 Adulyasak Y Cordeau JF Jans R et al 2015b The production routing problem A review of formulations and solution algorithms Computers Operations Research 55 141152 Ahuja R K 2017 Network flows theory algorithms and applications Pearson Education Ahuja R K Orlin J B Sharma D et al 1998 New neighborhood search structures for the capacitated minimum spanning tree problem Ahuja R K Orlin J B Sharma D et al 2000 Very largescale neighborhood search International Transactions in Operational Research 745 301317 Aiex R M Resende M G Ribeiro C C 2002 Probability distribution of solution time in GRASP An experimental investigation Journal of Heuristics 83 343373 Aiex R M Resende M G Ribeiro C C 2007 TTT plots a perl program to create timetotarget plots Optimization Letters 14 355366 Archetti C Bertazzi L Laporte G Speranza M G 2007 A branchandcut algorithm for a vendormanaged inventoryrouting problem Transportation Science 413 382391 Archetti C Bertazzi L Paletta G Speranza M G et al 2011 Analysis of the maximum level policy in a productiondistribution system Computers Operations Research 3812 17311746 Archetti C Speranza M G 2016 The inventory routing problem the value of integration International Transactions in Operational Research 233 393407 Armentano V A Shiguemoto A L Løkketangen A et al 2011 Tabu search with path relinking for an integrated productiondistribution problem Computers Operations Research 388 11991209 Avci M Yildiz S T 2019 A matheuristic solution approach for the production routing problem with visit spacing policy European Journal of Operational Research 2792 572588 Avci M Yildiz S T 2020 A mathematical programmingbased heuristic for the production routing problem with transshipments Computers Operations Research 123 Article 105042 Bard J F Nananukul N 2010 A branchandprice algorithm for an integrated production and inventory routing problem Computers Operations Research 3712 22022217 Bayley T Bookbinder J H 2019 A branchandcut algorithm for the production routing problem with coordinated and capacitated lotsizes and backorders In International workshop on lotsizingiwls2019 p 26 BeloFilho M Amorim P AlmadaLobo B et al 2015 An adaptive large neighbourhood search for the operational integrated production and distribution problem of perishable products International Journal of Productions Research 5320 60406058 Benders J F 1962 Partitioning procedures for solving mixedvariables programming problems Numerische Mathematik 41 238252 Bianchessi N Mansini R Speranza M G et al 2018 A branchandcut algo rithm for the Team Orienteering Problem International Transactions in Operational Research 252 627635 Bitran G R Tirupati D 1993 Hierarchical production planning Handbooks in Operations Research and Management Science 4 523568 Boudia M Louly M A O Prins C et al 2007 A reactive GRASP and path relinking for a combined productiondistribution problem Computers Operations Research 3411 34023419 Brahimi N Aouam T 2016 Multiitem production routing problem with back ordering a MILP approach International Journal of Productions Research 544 10761093 Brown G Keegan J Vigus B Wood K 2001 The kellogg company optimizes production inventory and distribution Interfaces 316 115 Caserta M Voß S 2009 Metaheuristics intelligent problem solving In Matheuristics pp 138 Springer Çetinkaya S Üster H Easwaran G Keskin B B 2009 An integrated out bound logistics model for FritoLay Coordinating aggregatelevel production and distribution decisions Interfaces 395 460475 Chandra P 1993 A dynamic distribution model with warehouse and customer replenishment requirements Journal of the Operational Research Society 681692 Computers Industrial Engineering 171 2022 108369 15 A Fortes et al Chitsaz M Cordeau JF Jans R 2019 A unified decomposition matheuristic for assembly production and inventory routing INFORMS Journal on Computing 311 134152 Clarke G u Wright J W 1964 Scheduling of vehicles from a central depot to a number of delivery points Operations Research 124 568581 Coelho L C Cordeau JF Laporte G et al 2013 Thirty years of inventory routing Transportation Science 481 119 Coelho L C Laporte G 2013 A branchandcut algorithm for the multi product multivehicle inventoryrouting problem International Journal of Productions Research 512324 71567169 Côté J Guastaroba G Speranza M 2017 The value of integrating loading and routing European Journal of Operational Research ISSN 03772217 2571 89105 Darvish M Archetti C Coelho L C 2018 Tradeoffs between environmental and economic performance in production and inventoryrouting problems International Journal of Production Economics Dayarian I Desaulniers G 2019 A branchpriceandcut algorithm for a productionrouting problem with shortlifespan products Transportation Science 533 829849 DíazMadroñero M Peidro D Mula J et al 2015 A review of tactical optimiza tion models for integrated production and transport routing planning decisions Computers Industrial Engineering 88 518535 Dolan E D Moré J J 2002 Benchmarking optimization software with performance profiles Mathematical Programming 912 201213 Fang X Du Y Qiu Y et al 2017 Reducing carbon emissions in a closedloop production routing problem with simultaneous pickups and deliveries under carbon capandtrade Sustainability 912 2198 Golsefidi A H Jokar M R A 2020 A robust optimization approach for the productioninventoryrouting problem with simultaneous pickup and delivery Computers Industrial Engineering Article 106388 Groër C Golden B Wasil E et al 2010 A library of local search heuristics for the vehicle routing problem Mathematical Programming Computation 22 79101 Hansen P Mladenović N 2001 Variable neighborhood search Principles and applications European Journal of Operational Research 1303 449467 Hansen P Mladenović N Todosijević R Hanafi S 2017 Variable neighborhood search basics and variants EURO Journal on Computational Optimization 53 423454 Kytöjoki J Nuortio T Bräysy O Gendreau M 2007 An efficient variable neighborhood search heuristic for very large scale vehicle routing problems Computers Operations Research 349 27432757 Li Y Chu F Chu C Zhu Z 2019 An efficient threelevel heuristic for the largescaled multiproduct production routing problem with outsourcing European Journal of Operational Research 2723 914927 Li Y Chu F Côté JF Coelho L C Chu C 2020 The multiplant perishable food production routing with packaging consideration International Journal of Production Economics 221 Article 107472 LópezIbáñez M DuboisLacoste J Cáceres L P Birattari M Stützle T 2016 The irace package Iterated racing for automatic algorithm configuration Operations Research Perspectives 3 4358 Lourenço H R Martin O C Stützle T 2003 Iterated local search In Handbook of metaheuristics pp 320353 Springer Lourenço H R Martin O C Stützle T 2019 Iterated local search Framework and applications In Handbook of metaheuristics pp 129168 Springer Lysgaard J 1997 Clarke wrights savings algorithm Department of Management Science and Logistics the Aarhus School of Business Miranda P L Morabito R Ferreira D et al 2018 Optimization model for a production inventory distribution and routing problem in small furniture companies TOP 261 3067 NevesMoreira F AlmadaLobo B Cordeau JF Guimarães L Jans R 2019 Solving a large multiproduct productionrouting problem with delivery time windows Omega 86 154172 Penna P H V Subramanian A Ochi L S 2013 An iterated local search heuristic for the heterogeneous fleet vehicle routing problem Journal of Heuristics 192 201232 Pisinger D Ropke S 2010 Large neighborhood search In M Gendreau J Y Potvin Eds International Series in Operations Research Management Science Handbook of metaheuristics 146 pp 399420 Springer New York Dordrecht Heidelberg London Springer Pochet Y Wolsey L A 2006 Production planning by mixed integer programming Springer Science Business Media Qiu Y Ni M Wang L Li Q Fang X Pardalos P M et al 2018 Produc tion routing problems with reverse logistics and remanufacturing Transportation Research Part E Logistics and Transportation Review 111 87100 Qiu Y Qiao J Pardalos P M et al 2017 A branchandprice algorithm for production routing problems with carbon capandtrade Omega 68 4961 Qiu Y Wang L Fang X Pardalos P M Goldengorin B et al 2018 Formulations and branchandcut algorithms for production routing problems with time windows Transportmetrica A Transport Science 122 Qiu Y Wang L Xu X Fang X Pardalos P M et al 2018 A variable neighborhood search heuristic algorithm for production routing problems Applied Soft Computing 66 311318 Reimann M Tavares Neto R Bogendorfer E et al 2014 Joint optimization of production planning and vehicle routing problems A review of existing strategies Pesquisa Operacional 342 189214 Reyes A Ribeiro C C 2018 Extending timetotarget plots to multiple instances International Transactions in Operational Research 255 15151536 Russell R A 2017 Mathematical programming heuristics for the production routing problem International Journal of Production Economics 193 4049 Schenekemberg C M Scarpin C T Pecora Jr J E Guimarães T A Coelho L C 2021 The twoechelon productionrouting problem European Journal of Operational Research 2882 436449 Shaw P 1998 Using constraint programming and local search methods to solve vehicle routing problems In International conference on principles and practice of constraint programming pp 417431 Springer Solyalı O Süral H 2017 A multiphase heuristic for the production routing problem Computers Operations Research 87 114124 Souza M J Coelho I M Ribas S Santos H G Merschmann L H d C 2010 A hybrid heuristic algorithm for the openpitmining operational planning problem European Journal of Operational Research 2072 10411051 Stålhane M Andersson H Christiansen M Fagerholt K 2014 Vendor managed inventory in tramp shipping Omega ISSN 03050483 47 6072 Thompson P M Orlin J B 1989 The theory of cyclic transfers In OPERATIONS RESEARCH CENTER WORKING PAPER MIT Toth P Vigo D 2014 Vehicle routing problems methods and applications SIAM Watanabe H Cirino R Soler W Santos M et al 2017 Solution methods for the integrated production routing problem Salesian Journal on Information Systems 119 5361