The cyclic production routing problem

This paper introduces the Cyclic Production Routing Problem (CPRP). The CPRP is an extension of the well-known NP-hard Production Routing Problem (PRP), which is a hard-to-solve combinatorial optimisation problem with numerous practical applications in the field of freight transportation, logistics and supply chain management. Under the PRP setting, a manufacturer is responsible for determining production decisions, as well as the timing and quantity of replenishment services offered to a set of geographically dispersed customers over a multi-period time horizon. The problem calls for jointly optimising the production, inventory, distribution and routing decisions. In this paper, the basic PRP model is modified to generate repeatable cyclic production and delivery schedules. A two-commodity flow formulation is proposed along with valid inequalities. Extensive comparisons between the basic PRP and the proposed cyclic variant on well-known benchmark instances are provided. The new variant is significantly harder to solve, especially when the vehicle fleet is limited. From a managerial perspective, the generation of cyclic production-routing schedules significantly increases all costs, whereas the number of vehicle routes required to implement a cyclic schedule is higher.


Introduction
The globalisation of markets and the increasing interdependence of economies underline the importance of efficient supply chain management.A supply chain is a network of organisations, resources and activities that are required for manufacturing, distributing and transporting a product.It includes all activities such as the production, inventory control, distribution and routing that may take place in geographically different areas.These activities are interdependent and determine the costs of the whole procedure.Therefore, their successful coordination promotes cost savings and increased effectiveness, contributing to the supply chain sustainability.
Despite the evident importance of jointly optimising all supply chain processes, the common practice is to tackle them individually in pursuit of simplicity and applicability.In particular, the most common optimisation approach is hierarchical starting from the production to consumption, without information being passed up the chain.Practically, this means that the master optimisation problem is divided into isolated sub-problems and the decision-making process becomes myopic, inelastic and often ineffective.For example, in a hierarchical optimisation approach of a productiondistribution system, the lot-sizing production decisions are made first, to minimise production and inventory costs.Then, the distribution and routing decisions are made with respect to the generated production plan, and no alterations of the production level decisions is allowed.A recent review concludes that integrating production, inventory and routing decision under a unified framework, leads to cost savings of 11.08% (Hrabec, Hvattum, and Hoff 2022).Regardless of the lower quality of this fragmented decision-making approach, it remains the common practice, because integrated decision making often leads to extremely complex logistics management problems that are difficult to frame and tackle (Archetti and Speranza 2016).
The scientific advances of the last years have enabled researchers and practitioners to formulate and solve more realistic problems that better reflect the business environment requirements.One of the most characteristic integrated problems arising in numerous supply chains, is the well-known Production routing Problem (PRP).The PRP is a hard-to-solve NP-hard combinatorial optimisation problem that calls for the joint optimisation of production, inventory, and routing decisions over a specific time horizon.A supplier is responsible for replenishing the inventories of geographically dispersed customers ensuring that no stock-outs occur at any period of the considered time horizon.Both customers and the depot/production facility face a per period inventory unit holding cost and have a maximum storage capacity.The problem incorporates the naturally cost-saving Vendor Managed Inventory (VMI) policy according to which the supplier is responsible for determining the quantities and the timings of customer replenishment visits, as well as the routes serving the customers.In that sense, PRP generalises the Inventory routing Problem (IRP) and incorporates the well-known lot-sizing problem (LSP) and the vehicle routing problem (VRP).The adoption of VMI policy creates a win-win situation, as on one hand, customers save resources from inventory control and order placement and on the other hand, suppliers are able to optimise production, inventory and distribution processes (Archetti and Speranza 2016).
The combined nature of the PRP is suitable for various industry fields and different applications.In practice, companies such as Kellogg (Brown et al. 2001) and Frito-Lay (Sila et al. 2009) have recorded significant savings after jointly optimising production and distribution operations.Inspired by the requirements and the particularities of different industries and sectors, researchers have introduced and studied several PRP variants.Motivated by a food company that distributes fresh meat to a network of stores in China, Qiu, Qiao, and Pardalos (2019) consider the PRP for perishable products and experiment with different selling policies to minimise the value losses.Also for the case of perishable products, Chan et al. (2020) extend the PRP in the context of sustainable food supply chain.They consider four objectives: the classic minimisation of the total system costs, as well as, the maximisation of the average food quality, the minimisation of the CO 2 emissions in transportation and production and finally, the minimisation of the total weighted delivery lead time.The perishability feature is also considered by Ghasemkhani et al. (2021).The authors formulate a multi-perishable product and multiperiod PRP with heterogeneous fleet, time windows and fuzzy parameters.The objective of the proposed model is to maximise the total profit (selling revenue reduced by the aggregation of the holding, production, transportation and utility preference costs).Motivated by a real case, Dayarian and Desaulniers (2019) model and solve the PRP of a catering company in Montreal that delivers meals with short life-span.The authors take into account business requirements such as multitrips and time-windows.With respect to the petrochemical industry, Schenekemberg et al. (2021) study the Two-Echelon PRP with pickups and inventory control of ethanol from the suppliers, production and inventory control of pure and commercial gasoline, as well as deliveries of commercial gasoline to the final customers in South American countries.Most recently, Farghadani-Chaharsooghi et al. (2021) integrate the PRP model with workforce planning for the case of a company that processes and delivers organic fruits and vegetables.They consider stochastic times, perishable products, and quantify the effects of workforce planning on costs and productivity offering interesting managerial insights.The interested reader is referred to the surveys of Díaz-Madroñero, Peidro, and Mula (2015) and Adulyasak, Cordeau, and Jans (2015a) on optimisation models which integrate production and transportation decisions.
An important aspect of real-world operational problems is the ability of an optimised set of decisions to be periodically applied in the long run, ensuring reduced variability and satisfactory performance (Grzegorz et al. 2021).It has been shown that for specific cases, the cyclic production planning can achieve a specified service target at lower total costs than a more costefficient but not non-cyclical approach (Nyen et al. 2009).Despite the applicability of several PRP model variants on real-world cases, to the best of our knowledge, the feature of cyclicity has not been previously studied.In this context, the term cyclicity is used to describe the capability of repeating the production, distribution and routing plans for several cycles.An optimal solution of the basic PRP model (Adulyasak, Cordeau, and Jans 2014a) assumes zero final inventory for both the production facility and customers (according to the basic PRP, a nonzero inventory would indicate sub-optimal solution due to the incurred inventory holding costs).Therefore, at the end of the considered time horizon, the whole supply chain network (production site and customers) is out of products.Thus, the generated optimised schedules cannot be repeated in a cyclic manner: in fact, stock-outs are inevitable right after the end of the examined time horizon.To overcome this limitation, the classic PRP model can be revisited through the prism of a cyclic time horizon to ensure that the final inventory levels of the depot and the customers are equal to the starting ones.In that sense, the production and distribution schedules along with the routing plans can be repeated for several cycles of length equal to the considered time horizon provided that the various problem parameters remain unchanged.Recently, Manousakis (2021) proposed a formulation for a periodic variant of the PRP that considers fixed equal starting and final inventories.
The present paper introduces the Cyclic Production routing Problem (CPRP), a practical PRP variant that calls for the optimal production, inventory and distribution plans which can be periodically applied in cycles.To the best of our knowledge, this is the first time that a cyclicity feature is integrated into a PRP model.The proposed model builds on and extends the work of Manousakis (2021) which considers a less flexible periodic PRP variant.It is formulated as a mixed integer linear programming (MILP) problem, strengthened with new and adapted from the PRP literature valid inequalities.Thanks to a flow formulation which omits the vehicle index, the MILP can be directly solved via a Branch-and-Cut algorithm (BnC) by any mathematical programming solver.We perform extensive computational experiments to compare the solutions obtained by the basic PRP and the proposed CPRP model on well-known benchmark instances.These experiments demonstrate that the more realistic CPRP is much harder to solve, especially when limited vehicle fleet is considered.Tested on benchmark instances, the requirement for cyclicity which leads to the production and distribution of extra product quantities significantly increases the total costs from 33% up to 56%, depending on the instance size.Therefore, in practice PRP may become misleading for a decisionmaker as it significantly underestimates the actual costs.Further comparative results are also provided regarding the required number of vehicles, as well as the average inventory levels of the depot and the customers.
The remainder of the paper is organised as follows.Section 2 provides the PRP related literature.In Section 3, the two-commodity flow CPRP formulation is given, together with adapted and new valid inequalities.Next, Section 4 reports extensive computational experiments and discusses managerial implications of incorporating the requirement for cyclicity into the PRP model.Finally, conclusions are summarised in Section 5.

Related literature
In this section, we discuss several studies which are relevant to our research.Initially, classic periodic vehicle routing problems are presented.Following this, the more flexible and hard-to-solve periodic and cyclic inventory routing problems are discussed.Finally, the section is concluded by reviewing literature related to the production routing problems.

Periodic vehicle routing problems
The well-known Vehicle routing Problem (VRP) (Toth and Vigo 2014) may be considered the cornerstone of optimisation for logistics operations.Vehicle routing problems seek to optimally assign customer orders to vehicles and sequence the service visits performed by each vehicle.Typically, the goal is to minimise the total distance travelled by the vehicles, while ensuring that all customer orders are served.Basic vehicle routing problems involve a single period.In practice, it is important to also consider different planning horizons for transportation activities: periodic deliveries, i.e. deliveries that are repeated in cycles of equal length are found in numerous fields and applications (Francis, Smilowitz, and Tzur 2008).On this basis, the Periodic Vehicle routing Problem (PVRP) Christofides and Beasley (1984) extends the classic VRP by considering a multi-period horizon during which the customer visits take place.Campbell and Wilson (2014) define the basic PVRP as the problem of selecting one of the given feasible visiting schedules of each customer (since each customer is not visited in every period), such that the total transportation cost is minimised.According to the same work, the delivery options for each customer found in the PVRP literature may be classified to three categories: (i) predetermined set of candidate visit schedules, (ii) requirement for equally distanced customer visits and (iii) imposing minimum and maximum numbers of periods allowed between two consecutive customer visits.For all cases, given a customer visit schedule, the delivery quantity of every customer visit is fixed.
The assumption of predetermined and inflexible visit schedule options and fixed delivery quantities limit the applicability of the PVRP model to real-world problems.An interesting generalisation of the PVRP is the Flexible Periodic Vehicle routing Problem (FPVRP) (Archetti, Fernández, and Huerta-Muñoz 2017).Under the FVRP setting, each customer has a known total demand that must be met through the planning horizon.A limit is imposed on the maximum quantity that can be delivered on each customer visit.As opposed to the more constrained PVRP, the FPVRP does not consider fixed service frequencies or fixed delivery quantities, and therefore it provides extra room for cost-savings.Regarding handling customer deliveries, FPVRP is similar to the well-known Inventory routing Problem (IRP) introduced in Bell et al. (1983): both decisions on when and how much to deliver have to be made (Mor and Speranza 2020).

Inventory routing problems
The family of IRPs is a broad class of typically multiperiod problems with numerous applications in several sectors and industries.The IRP calls for jointly determining the timing and quantity of customer deliveries, as well as the minimum cost vehicle routing in pursuit of the optimal coordination of inventory holding and transportation activities.By definition, IRP considers the VMI inventory replenishment policy, i.e. the supplier is responsible for delivering quantities to the customers, so that no stock-outs occur.The VMI policy commonly substitutes the short-sighted Retailer Managed Inventory policy under which each customer is responsible for placing orders.The VMI policy is one of the most important supply chain strategies followed by 84% of the companies with over a billion revenue (van den Bogaert and Jaarsveld 2021).VMI allows the vendor to effectively coordinate delivery activities, in pursuit of distribution cost savings.On the other hand, customers receive cost incentives and save time and effort on inventory management.According to Archetti and Speranza (2016), the inventory and routing cost savings achieved by the VMI policy may range up to approximately 10% for well-known benchmark data sets.VMI systems were introduced in the literature in the context of liquified air products distribution.Since then, several VMI systems have been implemented for several sectors: automotive industries, electronics assembly, chemicals industries, vending machines for juice or foods, chain stores, maritime logistics and many other (Andersson et al. 2010).The interested reader is referred to the survey of Coelho, Cordeau, and Laporte (2014) for an insightful look on the IRP literature.
The role of periodicity in the sense of repeated distribution patterns has also been examined under the IRP setting.In comparison to the classic IRP model (Archetti et al. 2014), the Periodic Inventory routing Problem (PIRP) incorporates additional constraints that enforce equality of the initial inventory levels and the inventory levels at the end of the planning horizon for both the depot and the customers.Therefore, the generated distribution schedule may be repeated for as long as the input parameters remain unchanged.This approach enables solving problems without optimising over an infinite horizon (van Anholt et al. 2016).The periodicity feature strengthens the practical applicability of PIRP: For example, Gaur and Fisher (2004) solve the problem of scheduling and routing replenishment operations for a supermarket network.The optimised schedule is repeated weekly.Over the first year of the schedule implementation, the distribution cost savings were equal to 4%.The Selective and Periodic IRP (SPIRP) is solved by Aksen et al. (2014) via an adaptive large neighbourhood search algorithm.It calls for the generation of a distribution plan which is repeated on a weekly basis.The examined model is faced by a company which collects used vegetable oil from sources, such as restaurants and hotels and reuses it to produce biodiesel.In contrast to the basic IRP version under which all customers must be visited to prevent stock-outs, this variant allows selecting which customers to visit based on the profitability.Similarly, Montagné, Gamache, and Gendreau (2019) develop a constructive algorithm based on shortest path and split procedures for a real case of reusable waste oil collection in Canada.Interestingly, the tests performed on real-world problem instances of up to 3,000 customers served on a 30-day time horizon, show that the algorithm manages to outmatch the cost-effectiveness of the actual company solution by up to 20%.
In the case of PIRP, the requirement for repeatable replenishment schedules is met by enforcing same inventory levels at the start and the end of the studied time horizon.Therefore, the whole delivery schedule may be repeated in a cyclic manner, with the cycle time being equal to the time horizon considered (e.g., a week).Other IRP approaches, incorporate the periodicity feature by allowing different cycle times for each customer or for each route.For instance, a customer with high demand and high proximity to the depot may be served every second day, whereas, a customer with low demand rate and far from the depot may be served once per week.Similarly, a route may be repeated every day and another route every third day.In these cases, different delivery patterns are used for each customer and the whole schedule may be repeated with a cycle time equal to the lowest common multiple of the customer, or the route cycle times.On this basis, another IRP variant which is referred to as the Cyclic IRP (CIRP) is introduced.The objective of CIRP is to find a cyclic distribution pattern that minimises the long term transportation and inventory costs.The cycle time is an important decision variable and it is often dictated by economic order quantity (EOQ) models.Aghezzaf, Raa, and Van Landeghem (2006) propose a long-term IRP model with constant demands and consider economic order quantity (EOQ) policies for inventory management.A column generation based approximation method for solving the non-linear formulation is proposed.Raa and Aghezzaf (2009) develop a practical solution approach that considers cargo handling times, customer inventory capacities, and maximum driving limits for vehicles.Chitsaz, Divsalar, and Vansteenwegen (2016) propose a two-phase iterative procedure which consists of two heuristic methods: the first one produces routes, whereas the second one combines and schedules routes to generate cyclic distribution plans.A two phase approach is also proposed by Raa and Dullaert (2017): in the first phase, cycle times are chosen and the routes are designed, whereas in the second phase, the routes are assigned to vehicles and the cycle times are adjusted, to minimise the number of vehicles used.The CIRP model is applicable to various industries.For example, Zenker, Emde, and Boysen (2016) study the problem of producing cyclic tours when customers are located along a line.This operational scenario occurs in liner shipping (feeder ships that service inland ports along a stream) and facility logistics (tow trains which deliver bins to the stations of an assembly line).More recently, Bertazzi et al. (2020) study a problem where components are collected from a set of supplier locations and delivered to a manufacturing plant in cycles.The authors present polyhedral studies of the convex hull of the problem and propose a branch-and-cut algorithm to solve the provided problem formulation.Interestingly, computational experiments show that the cyclic formulation is significantly harder to solve compared to the standard non-cyclic IRP formulation.More specifically, optimal solutions are found for problem instances of up to 25 customers with the cyclic model, whereas for the non-cyclic version optimal solutions are found for problems of up to 50 customers when a single vehicle is considered.

Production routing problems
The Production routing Problem introduced by Chandra (1993) extends and generalises the IRP by considering production decisions.It has been pointed out that an integrated production, inventory and routing planning can reduce the total operating cost by 3-20% (Chandra and Fisher 1994).As a result, PRP has received increasing research attention, with most works published during the last few years.
Due to the high complexity of the PRP, the majority of early approaches are based on decomposing the master problem to subproblems that are usually tackled sequentially via metaheuristics such as Greedy Randomized Adaptive Search Procedure (GRASP) with path relinking (Boudia, Louly, and Prins 2007), two-stage local search algorithms (Boudia, Aly Ould Louly, and Prins 2008), memetic algorithms with population management (Boudia and Prins 2009), or tabu search schemes with path relinking and both short and long-term memories (Armentano, Shiguemoto, and Løkketangen 2011).The exact approaches are scarce and capable of tackling relatively small problem instances.For example, for the single vehicle case, Archetti et al. (2011) apply an exact approach to solve instances of six time periods, whereas Qiu et al. (2018b) propose a Branch-and-Cut (BnC) method with new valid inequalities which can solve instances of only three periods and up to 100 customers.Adulyasak, Cordeau, and Jans (2014a) compare the performance of vehicle-indexed and nonvehicle indexed formulations for several test problems.To do so, single core and parallel computing BnC algorithms are developed and applied.The same relatively small instances are solved for perishable inventory by a BnC method strengthened with lot-sizing and lifted Miller-Tucker-Zemlin subtour elimination constraints (Qiu, Qiao, and Pardalos 2019).Adulyasak, Cordeau, and Jans (2015a) introduce a Benders decomposition approach which is implemented in a single branch-andbound tree and enhanced through lower-bound lifting inequalities, scenario group cuts, and Pareto-optimal cuts for a stochastic IRP variant with demand uncertainty.Brahimi and Aouam (2016) present a mixed integer linear programming formulation for the production routing problem with backordering.More recently, Zhang et al. (2021) propose a Benders decomposition approach, to solve the PRP with order-up-to-level (OU) replenishment policy (every time a customer is visited, the quantity delivered is such that the maximum inventory level is reached).
The most common and efficient approaches for hardto-solve routing problems combine metaheuristics and mathematical programming methods and are referred to as matheuristics (Archetti and Speranza 2014).Matheuristics have shown merit for several PRP variants.Adulyasak, Cordeau, and Jans (2014b) introduce an optimisation based adaptive large neighbourhood search heuristic: the binary production and setup decisions are handled by an enumeration scheme and upperlevel search operators, whereas all continuous decisions related to production, inventory and delivery decisions are made by solving a minimum cost network flow subproblem.Decomposition approaches are also proposed by Absi et al. (2015) and Miranda et al. (2018): the authors employ two-phase iterative schemes with the production-distribution and routing subproblems solved by Traveling Salesman Problem (TSP) and Vehicle routing Problem (VRP) solvers, respectively, within a local search framework.In the field of genetic algorithms, Senoussi et al. (2018) use decomposition techniques and solve MIPs to generate good quality individuals from the population.
The basic PRP version (Adulyasak, Cordeau, and Jans 2015a) has served for benchmarking several highquality solution approaches, most of which are matheuristic algorithms: Russell (2017) uses mathematical programming for a relaxed PRP version, to determine an initial solution.This solution is completed by a tabu search scheme based on the concept of seed routes.Solyalı and Süral (2017) propose a five-phase heuristic, with overlapping subproblems which are formulated as MIPs and solved via exact algorithms.In a similar manner, other multi-phase approaches decompose the problem into the setup, distribution and routing decision levels which may be tackled either by exact algorithms (Chitsaz, François Cordeau, and Jans 2019), or with fix-and-optimise strategies (Li et al. 2019).On the same basis, Avci and Yildiz (2019) propose a decomposition of the main problem into subproblems to reduce its complexity.The distribution and routing subproblems are solved heuristically (iterated local search), whereas the lot-sizing, inventory and delivery quantity decisions are handled by solving the corresponding MIP models.Hybrids of metaheuristic algorithms have been also introduced.Qiu et al. (2018a) adopt a skewed general variable neighbourhood search algorithm for the delivery schedule and a guided variable neighbourhood descent algorithm for the routing subproblem, whereas the production and inventory quantities are obtained by solving a production-inventory MIP subproblem.Most recently, Schenekemberg et al. (2021) introduced a parallelised hybrid algorithm that combines local search procedures within a traditional BnC scheme.

Mathematical formulation
In this section, the mathematical formulation of the proposed Cyclic Production routing Problem is provided.Several valid inequalities are adopted from the PRP literature.

The cyclic production routing model
The Cyclic Production routing Problem (CPRP) seeks production, distribution, inventory and routing plans which can be repeated in a cycles of length equal to the studied time horizon.As per the basic PRP model, a product supplier is responsible for replenishing the inventories of customers over a fixed cyclic planning horizon which consists of a predetermined set of periods and guarantees that no stock-outs occur.The supplier decides the time and the quantity of production and distribution, as well as the routes serving the customers.The differentiating factor of the introduced CPRP variant is the requirement of equal ending and starting inventories for both the production facility and the customers.As the model is clearly affected by the initial inventory levels, the starting inventories of both the customers and the production facility are considered to be decision variables.Unlike, the work of Manousakis (2021) for a periodic PRP which considers starting inventories as parameters, handling initial inventories as variables offers flexibility and increased cost savings.The cost of preparing and establishing the optimal inventory levels is considered negligible.If this is not the case, this preparation cost may be reflected in the objective function.The twocommodity flow PRP model of Manousakis et al. (2022) is adapted to the examined CPRP.It is worth mentioning that this type of modelling has shown promising results for the IRP which can be regarded as a special case of the basic PRP (Manousakis et al. 2021).
Let G = (V, E) be an undirected graph where V = {0, 1, . . ., n, n + 1} is the set of nodes and E = {(i, j) : i, j ∈ V, i < j} is set of edges.Node 0 represents the production facility and node n + 1 represents a production facility clone, where the inverse flow originates from.The set of customers is denoted as V C = {1, 2, . . ., n}.The considered planning horizon is considered to repeat in cycles.Each cycle consists of |T| periods (so that the plans are repeated every |T| periods).Let T = {1, 2, . . ., |T|} denote a complete cycle of the planning horizon.To facilitate the model description, let prev(t) denote the index of the period before period t.This is prev For the first period of the cycle (t = 1), prev(1) = |T|.
Each node i ∈ V (both the production facilities and the customers) incurs a unit holding cost h i for every period t ∈ T. Customer i ∈ V C faces a per period t ∈ T demand d ti and has a limited maximum inventory capacity equal to L i .At the start of period t ∈ T the production facility may choose to produce any nonnegative quantity p t up to the production capacity limit C. At any period t ∈ T, if p t > 0, then a production set up cost s t is incurred.In addition, a per product unit production cost of u t also applies.The produced quantity may be used directly at the same period t to satisfy customer demands.Each edge (i, j) ∈ E is associated with a non-negative travel cost c ij that represents the cost of a vehicle for traversing this edge.The supplier delivers any non-negative quantity to every customer i ∈ V C at t ∈ T such that no stock-outs occur.For every period, a homogeneous vehicle fleet of |K| vehicles K = {1, 2, . . ., |K|} each of capacity Q is available at the depot.
We assume a symmetric transportation cost matrix, i.e. ∀i, j ∈ E, c ij = c ji , that satisfies the triangle inequality.Binary undirected routing variables x tij = 1 for i, j ∈ V, i < j take value 1, iff any vehicle k ∈ K traverses edge (i, j) in any direction at period t ∈ T. Binary variables z ti are equal to 1, iff i is visited by any vehicle at period t ∈ T. Similarly, binary variables y t are 1, iff production takes place at the production facility during period t.Nonnegative continuous flow variables f tij with i < j and f tji with i > j represent the load and the residual capacity of the vehicle travelling from i to j at time t ∈ T, respectively.The quantities produced at each period t are denoted as p t , whereas q ti is the non-negative product quantity delivered to customer i ∈ V C at time t ∈ T. The inventory at the end of period t for every node i ∈ V is captured by the continuous variable I ti .According to the cyclicity constraint, the end inventories of each node i ∈ V have to be equal to the starting ones, i.e.I |T|,i = I 0i .The objective of the CPRP is to minimise the overall fixed and variable production, transportation and inventory holding costs for both the supplier and the customers.
Below, the two-commodity flow formulation for the CPRP with the maximum level (ML) inventory policy is provided, henceforward referred to as CPRP: x ti,(n+1) t ∈ T (4) p t ≤ Cy t t ∈ T (9) q ti t ∈ T (10) The objective function (1) represents the sum of the depot (production facility) setup and variable production costs, the inventory holding costs of both the depot and the customers over the planning horizon, and the transportation costs.Note that no inventory holding costs are considered for the depot clone node.Constraints (2) are the degree constraints for the customers, whereas constraints (3) and ( 4) bound the number of vehicles leaving and returning to the depot, respectively.Constraints (5) ensure that the total flow (normal and inverse of each edge) is equal to the vehicle capacity, if and only if the edge is traversed by a vehicle.Constraints (6) are the flow continuity constraints that ensure that the total sum of the flows originating from a node are reduced by the delivery quantity absorbed by the node.If any customer i ∈ V C , is not visited at period t (z ti = 0), then constraints (2) and ( 5) make sure that both the routing and the flow variables to and from this node are zero.As a result, constraints (6) force the delivered quantity to be zero, q ti = 0. Thus, the delivered quantity to a customer may be positive, iff a visit is performed at this customer z ti = 0.The total product quantity starting from the depot is given by ( 7).Constraints (8) ensure that no products return to the production facility, or in other words make sure that the total product quantity leaving the depot is delivered to the customers.Constraints (9) allow production to take place, only when the production facility is open.In this case, the produced quantity cannot exceed the total production capacity limit.Constraints ( 10) and ( 11) are the inventory flow preservation constraints over the periods of the planning horizon for the depot and the customers, whereas constraints (12) dictate the ML policy.Note that constraints ( 10) and ( 11) jointly set the inventory of every node at the end of the horizon equal to the starting inventory, and thus ensure the repeatability of the decisions made for the next cycle of the planning horizon.Finally, constraints ( 13)-( 19) enforce integrality, as well as lower and upper bounds on the decision variables.The flow variables f tij are defined by ( 16) for both directions of each edge, because the existence of flows imposes direction to the model.It should be noted that a worst-case analysis shows that the ratio of the optimal solution objectives under the cyclic and non-cyclic versions of PRP ( Z(CPRP) Z(PRP) ) is unbounded.The related proof is provided in Appendix A.
The proposed flow formulation eliminates subtours by jointly considering constraints ( 5)-( 8).Therefore, there is no need to introduce 2 |V| additional constraints to implement the classic DFJ sub-tour elimination constraints (i.e. one for each subset of V), or new variables for the alternative MTZ sub-tour elimination constraints.Hence, the proposed formulation can be solved as-is via any off-the-shelf Mixed Integer Linear Programming (MILP) solver.In the case that the initial inventories cannot be set to the desirable levels, the model can be modified to incorporate fixed initial inventories.Assuming a given initial inventory I init,i for each node i ∈ V, the following constraints have to be added to the model: Constraints ( 20) force the ending inventory of every node to be equal to the desired initial inventory level.
Regarding the maximum allowed delivery quantity of customer i at period t, we distinguish between two cases that have appeared in the production routing problem literature.The common practice, is to allow the delivered quantity q ti to exceed the maximum capacity L i making sure that the excessive quantity is consumed during this period t, so that the inventory I ti at the end of the period does not exceed maximum capacity L i .However, other research works, adopt a stricter assumption forcing the delivery quantity q ti ≤ L i , ∀i ∈ V C , t ∈ T (Adulyasak, Cordeau, andJans 2014a, 2015b;Qiu et al. 2018a).In the present work, we consider the more frequent former case.However, the proposed model can be modified to capture the latter case by replacing constraints ( 12) and ( 17) with ( 21) and ( 22):

Valid inequalities
This section presents valid inequalities that can be used to tighten the proposed CPRP formulation.The valid inequalities are adopted from Manousakis et al. (2022).
For the sake of brevity, only the adapted valid inequalities are presented below, whereas the valid inequalities that are used as-is to tighten the CPRP formulation are described in Appendix B.
The following valid inequalities bound the maximum delivery quantities: Similarly to constraints (17), valid inequalities (23) bound the maximum delivery quantities with respect to the maximum inventory capacity and the vehicle capacity.The difference is that they are multiplied by the visit variables, and thus ensure that the delivery quantity of a customer may be positive, only if this customer is visited.Note that this is also guaranteed by a combination of constraints in the basic formulation (see model description), but preliminary experiments showed that adding these inequalities straightforwardly, contributes to the lower bound of the root node relaxation.The concept of sub-deliveries initially introduced by Desaulniers, Rakke, and Coelho (2016) is used to impose a lower bound on the visits per customer for each period.Specifically, Lefever et al. (2018) bound the visits prior to a specific period t with the minimum number of subdeliveries (according to the customer capacity) that are required to satisfy the demand of period t.For example, a customer with maximum capacity L i = 60 and demand d ti = 20, has to be visited at least once during periods t ∈ {3, 4, 5, 6} in order to satisfy the demand for period t = 6.
The calculation of set P − it is described in detail in Lefever (2018) and can be directly applied to the examined CPRP model.It is therefore omitted for the sake of brevity.

Computational results
This section presents computational results obtained on well-known benchmark instances for the proposed Cyclic PRP model (CPRP) and the classic non-cyclic PRP model.The benchmark data sets are firstly discussed along with the compared production routing settings.Then, a graphical representation of the solutions obtained under the different formulations is given, to visualise the cyclicity effect on the solution structure.After this, the CPRP and PRP formulations are evaluated with respect to solvability and computational times.Next, we present an analysis of the cyclicity effect on the various objective categories (i.e.setup, production, inventory, routing), as well as, their contribution to the total solution cost.Following this, we discuss the effect of cyclic schedules on the inventory levels and the routing decisions.

Benchmark data sets
Our experiments are conducted on two of the most widely used PRP benchmark data sets.The first one was introduced by Archetti et al. (2011), henceforward referred to as data set A, whereas the second one was introduced by Adulyasak, Cordeau, and Jans (2014a), henceforward referred to as data set B.
Data set A consists of A1, A2 and A3 instance sets of six time periods and 14, 50, and 100 customers, respectively.For each customer, the per period demand is constant, the initial customer inventory is non zero, whereas the production facility has zero initial inventory.Additionally, no limitations to the production quantity and inventory capacity are imposed on the production facility.Each of the three sets includes 480 problems: five instances for each of 96 instance types (96 × 5 = 480).A1 considers one vehicle, whereas A2 and A3 assume unlimited vehicles.The 96 instance types are categorised to four main classes.Class I (1-24) corresponds to the instances with standard production, inventory and transportation costs.Class II (25-48) and Class III (49-72) are characterised by high production (i.e. unit production costs are multiplied with 10) and transportation variable costs (i.e.coordinates multiplied by 5), respectively.Class IV (73-94) considers no inventory holding cost for customers.
Data set B was created by considering various configurations for a subset of data set A instances.Its instances feature from two up to four vehicles and from 10 up to 50 customers.In general, they are more diverse in terms of the problem characteristics.
Note that for the cyclic version of the PRP, the initial inventories of both the production facility and the customers are ignored as they are decision variables.

Production routing models
In the case of CPRP, the requirement for cyclic routing and production schedules promote: (i) the production of the additional quantities during the last periods and (ii) the delivery of these quantities to the customers just before the end of the time horizon to avoid inventory the holding costs of previous periods.In particular, when the per period production quantities and the available vehicles are unbounded (as in the case of sets A2 and A3), it is expected that the optimal CPRP solutions involve the production and distribution of large quantities at the last period of the planning horizon which is unrealistic in practice.Hence, we consider a CPRP setting with limited vehicles named CPRP-LV : it is the CPRP model with the vehicle fleet parameter |K| fixed at the maximum number of vehicles used across all periods of the PRP solution.This setting is only applicable for data sets A2 and A3, since data sets A1 and B consider a given number of vehicles for each instance.
In brief, we compare the following problem settings: (1) CPRP: the Cyclic PRP model as described in Section 3.1.(2) PRP-Adj: the basic PRP model of Manousakis et al. (2022).The initial inventories of the depot and the customers are set equal to the initial inventories of the solutions obtained for the CPRP model which treats them as decision variables.
(3) CPRP-LV: the CPRP model for data sets A2 and A3 with limited vehicles.The fleet size |K| is fixed at the maximum number of vehicles used across all periods of the PRP-Adj solution (max t∈T j∈V C x t0j ).Note that by considering limited vehicle availability, if the starting inventories of CPRP are used, some instances may be infeasible although they can be solved for other starting inventories.Therefore, unlike PRP-Adj, the starting inventories are decision variables and not fixed equal to the ones of the CPRP.

Solution method and tuning
All models were coded in C# and solved via branch-andcut (BnC) with Gurobi 9.0.2.solver.A single thread was used to enable comparisons.The computational experiments were carried out on an Intel Core i7-7700 3.60 GHz CPU with 16 gigabyte RAM x64 Windows 10 machine.
In terms of Gurobi settings, the MIPGap was reduced to 10 −8 to accommodate objectives with many digits and the heuristic were set to 0.25 to find solutions in large instances, whereas all other parameters were set to default.Branching priority was given to the y variables.Finally, a two-hour CPU time limit per instance was imposed for data set A, whereas a three-hour CPU time limit was imposed for data set B, to ensure that the nine-period instances are solved.

Comparison of production routing problem variants: visual analysis of solutions for different models
This section presents a visual analysis of the solutions produced under the CPRP, the PRP-Adj and the CPRP-LV settings, for the example instance 70 of Class III of set A2.A similar analysis is provided in Appendix C for the MVPRP_C30_P9_V4_I3 instance of data set B.

Route plots
Figure 1 illustrates the routes for the CPRP (left), PRP-Adj (middle) and CPRP-LV (right).Each sub-figure consists of six plots which depict the routes for periods t = 1 to t = 6.Each colour represents a specific vehicle.Blue dots represent customer locations, whereas the black square is the production facility where all routes have to begin and finish.Under the CPRP setting, 15 routes in total are required to serve all customers, four and five of which are performed at periods t = 2 and t = 3, respectively.No visits take place for three out of six periods.The cyclicity feature implies that the final inventories are equal to the starting ones.Since holding costs apply and given that no constraint is imposed on the vehicle number or the production capacity, the optimal policy can be summarised as follows: produce the quantities of the final inventories and serve the customers during the last period.Therefore, more vehicles (six) are needed to serve customers during the last period.On the contrary, for the PRP-Adj setting fewer vehicle trips are planned (10), whereas no deliveries are made during the last two periods.The fact that empty inventories are allowed for the end of the planning horizon leads to significantly lower routing costs compared to the ones required by a repeatable solution (CPRP).The same holds for the required vehicle availability.Finally, it is clear than the PRP-Adj solution is front-loaded with routes in comparison to the more applicable CPRP solution.
Under the CPRP-LV setting, the maximum number of vehicles used at a single period is four, as for the PRP optimal solution.This vehicle restriction drastically affects the solution structure: almost all four vehicles are utilised during four periods.The quantities transported are larger, and therefore routes appear shorter (each route consists of fewer customers in comparison to the two aforementioned models).This solution seems more balanced and realistic and can be implemented with a smaller vehicle fleet.Obviously, the routing decisions are greatly affected by the introduction of cyclicity into the decision-making.

Inventory heat maps
Figure 2 illustrates the per period depot and customer inventory levels for all three settings, as a heat map.Rows correspond to time periods and the columns to the nodes.The first row depicts the starting inventory, whereas the first column represents the depot inventory level.The value of each cell is the percentage of capacity utilisation (100% indicates that the maximum inventory is stored).The darker colour indicates higher inventory.
The heat maps of the CPRP and PRP-Adj models are similar for the first periods of the planning horizon.Notice that the starting inventories are identical (the optimal starting inventories for the CPRP are given as input for the PRP-Adj model.For both settings, the inventories are high between periods t = 2 and t = 4, when they start to gradually reduce.However, for CPRP the required inventories are replenished during the last period due to the unlimited routes that are allowed.For PRP-Adj, the last row visualises the unrealistic assumption of zero end inventories: All nodes hold no stock, so that stock-outs at the next period are unavoidable.Both inventory heat maps of the CPRP and PRP-Adj solutions depict the limited applicability of these two models with respect to the unlimited vehicle availability and the zero ending inventories, respectively.On the contrary, the rightmost CPRP-LV solution inventory heat map is significantly different.Customer inventories do not follow a common pattern.Instead, each customer seems to experience peaks in inventory when the respective deliveries take place.This is due to the limited vehicle availability which make large delivery quantities necessary to ensure feasibility.

Comparative analysis
This section summarises the performance of the three compared settings for both data sets.Table 2 provides results grouped by the instance class, whereas the Table 3 results are grouped by the number of periods, as well as the vehicle fleet size.
Six columns are provided for each setting.Columns S, O, I and NS refer to the number of solved, optimally solved, infeasible and not solved (CPU time limit reached with no solution) instances (note that all the instances of data set B are solved, and therefore columns I and NS are omitted).Column Gap is the upper-lower bound gap upon termination calculated as 100 × UB−LB UB , where   UB is the best objective score and LB the lower bound.
Finally, column T provides the computational time in seconds.
Concerning data set A, the extra quantities needed to be produced under the cyclic settings render 160 out of the 480 A1 instances infeasible due to the limited capacity of the single vehicle.Moreover, the cyclic formulation appears harder to solve: for the PRP-Adj, 262 instances are solved optimally, whereas only 16 can be optimally solved within the same CPU time, under the cyclic setting.All 50-and 100-customer instances are solved, under the CPRP and the PRP-Adj models.Interestingly, under the CPRP-LV setting, a different performance is observed: Firstly, all A1 solutions are identical to the A1 CPRP solutions, because the fleet is already restricted to one vehicle.However, for the 50customer A2 instances, 7/480 are found infeasible.For the harder A3 set, 16/480 instances are not solved in time, whereas four test problems are found infeasible due to the constrained vehicle fleet size.In all cases, the gaps are much higher than the gaps of CPRP.The PRP-Adj model demonstrates much lower gaps compared to both cyclic variants.Concluding, the introduction of cyclicity clearly makes the formulation more difficult to solve and increases the optimality gaps and the required computational times.The effect is more evident when limited vehicle fleet availability is considered.Interestingly, the larger gaps are observed for instances of Class III which feature increased transportation costs.
Similar findings are observed for the second data set.The CPRP model is much more difficult to solve: larger computational times are required and larger final gaps are observed.Under the PRP-Adj setting 28 more instances are solved to optimality, whereas the gap in nearly 0.4% lower.Finally, the number of vehicles has a greater impact on the ability to obtain high-quality solutions compared to the number of periods, which seems to have a less noticeable effect.

Effect of cyclicity on cost categories
To better understand the effect of the cyclicity constraints on the solution objective, we examine the four cost categories considered in the objective function.Firstly, Table 4 provides an overview of the percent contribution of each cost category to the total cost for both data sets.Figure 3 presents the ratios of PRP-Adj and CPRP-LV solution cost over the CPRP solution cost which is used as the baseline for the comparisons (i.e.cost PRP−Adj cost CPRP and cost CRPR−LV cost CPRP ).The results are aggregated by the sets.Each plot illustrates the cost ratios observed for the various objective components: production setup cost (SC), unit production cost (PC), routing cost (RC), total inventory holding costs for both the supplier and the customers (IC), and the total cost (TC).The horizontal black line indicates the 100% value, i.e. cost equal to the baseline cost of the CPRP solution.The corresponding figure for the second data set is provided in the Appendix (Figure C3).Note that due to the single vehicle restriction of instance set A1, the CPRP and CPRP-LV solutions are identical.Moreover, only instances for which a solution is generated within the CPU time limit are considered.
As shown in the Table 4, the objective values are dominated by the unit production costs, for both data sets.The requirement for repeatable solutions forces extra quantity to be produced.As a result, the share of the production costs is significantly higher for the cyclic variants.The various costs for the data set B instances are more balanced, i.e. the setup, routing and inventory costs have a stronger contribution to the total cost.Additionally, the contribution of the setup costs in the total objective is higher compared data set A. On the contrary, for the instances of data set A, the routing costs account for a significant part of the total costs, whereas the setup and the inventory costs are not as impactful.
Concerning, Figure 3, the PRP-Adj setting leads to considerably lower total costs (65.64% of the CPRP model costs) with the greatest cost differences observed for the inventory and unit production costs.The results for set A2 show a marginal increase in the total cost for the CPRP-LV .The impact of limited vehicles is mainly seen in the increase of inventory costs by 23.18%.Again, the costs of the PRP-Adj solutions are on average equal to 64% of the corresponding CPRP solution costs.All cost categories of the PRP-Adj solutions are significantly lower, with the largest difference observed between the two settings being the unit production cost, due to the fact that less total quantity needs to be produced (since zero end inventories are allowed).The findings of set A3 are similar except for the inventory cost difference: the CPRP-LV costs are 182.45% of the CPRP costs.Also, there are higher setup costs which indicate that due to the limited vehicles and less frequent visits, the production facility needs to be opened for extra periods.Of course, this depends on the relationship between the production and the inventory holding costs.For the majority of instances, the inventory costs make up for only a small portion of the total solution cost.
Similar results have been obtained for data set B: the total costs for PRP-Adj are on average 17.34% lower than the corresponding CPRP total costs.

Inventory levels
As shown in 4.6, the cyclicity constraints are related to higher inventory costs (Figure 3).Further exploring the inventory aspects, we present box-and-whisker plots (Figure 4) with the average inventory level kept on customers and the production site, for the three compared model settings.Each data point refers to the average inventory level of a single benchmark problem.This average inventory level is calculated as the average inventory level over all nodes and periods, whereas the average inventory level for a given node i ∈ V period pair t ∈ T is equal to the percent utilisation of the maximum capacity, i.e. 100 × I ti L i .The whiskers extend to a maximum of 1.5 times the interquartile range.Finally, Table 5 contains the average inventory level for each setting, as well as, the total product units kept in inventory throughout the planning horizon, i.e. the sum of the per period inventories (i.e. a product unit that remains in stock for three days contributes with three units to the inventory sum.The corresponding figure and table for data set B are available in the Appendix (Figure C4 and Table C1).
The average inventory level and the sum of units kept in stock under the CPRP-LV is higher than those of the PRP-Adj.In addition, larger average inventory variance is observed among the benchmark instances for both cyclic settings compared to the PRP-Adj model.

Vehicle fleet and routes
In this section we examine the impact of the cyclicity requirement on the routing plans.Figure 5 contains a bar plot for each instance set.Each plot illustrates the number of routes for periods t = 1 up to t = 6, the total number of routes during the planning horizon, and the maximum (over all periods) number of routes for all three compared problem settings.The presented values are averages over all benchmark instances that were solved.The corresponding figure for the instances of data set B is given in the Appendix (Figure C5).For the A1 instances, the vehicle fleet is limited to one vehicle.The CPRP and CPRP-LV models yield identical solutions, whereas for the PRP-Adj model, a significant decrease on the average number of routes is observed for the last two periods.In addition, on average, the number of total vehicle trips for serving customers is lower compared to the cyclic versions.The differences between the three settings are clearer for the A2 and A3 sets: the total routes are approximately 15.70 ( A2) and 30.88 ( A3) for the PRP-Adj solutions, whereas for CPRP they increase to 21.15 ( A2) and 41.03 ( A3), respectively.Similarly, under the CPRP-LV setting, the total routes are 20.46 and 40.00, respectively.Regarding the distribution of the number of routes across periods, we observe that the CPRP has more routes than CPRP-LV on most periods.As expected, the PRP-Adj model requires significantly less routes, especially during the last periods of the planning horizon.Overall, the constrained vehicle fleet of the CPRP-LV setting forces routes to be distributed evenly among periods.From a different perspective, routes cannot be performed at the optimal timing with respect to the inventory holding cost, thus inventory needs to be accumulated at customer locations.
Figure 5 demonstrates the effect of considering a fixed vehicle fleet on the solution structure.In specific, when the vehicle fleet is constrained, deliveries take place consistently throughout the planning horizon, to meet the total customer demand.As a result, products tend to be accumulated at the customers (higher inventory levels) in pursuit of feasible solutions.On the contrary, when unlimited vehicles are assumed, the per period number of routes has greater variance through the planning horizon.For a practitioner, this implies that in order to implement the route schedules, there must be increased vehicle availability for meeting the requirements of the peak periods.This can be achieved either by a proprietary vehicle fleet or by outsourcing the delivery operations to a third party.If a proprietary fleet is assumed, then high idle times and maintenance costs, as well as low vehicle utilisation are unavoidable.On the other hand, if the increased number of vehicles is satisfied via outsourcing on specific periods, then the transportation costs are increased.Such realistic operational constraints have been studied in the context of routing problems: the concept of consistency has been introduced, to reflect some common standards for service quality.For instance, for the IRP, Coelho, Cordeau, and Laporte (2012) discuss six consistency features, including the vehicle filling rate: a minimum vehicle capacity utilisation is imposed.This is a common restriction used by enterprises, to avoid idle time costs and to balance driver workloads towards fairness.The authors find that the aforementioned minimum capacity utilisation constraint seems to be the most expensive consistency feature out of the six tested: the costs are increased by 5.88% up to 39.44% for the test problems considered.Similarly, our results indicate that incoporating realistic features such as cyclicity and limited fleet availability into the decision-making increases the costs for the production routing problems.It is worth mentioning that outsourcing has been studied for a multiple product PRP variant by introducing a per unit outsourcing cost (Li et al. 2019).
The data sets A and B are available by Archetti et al. (2011) and Adulyasak, Cordeau, and Jans (2014a), respectively, whereas the solution scores for the three settings are provided in the supplementary material.

Conclusions
This paper introduces a realistic variant of the wellknown production routing problem (PRP).Under the classic PRP, the network of customers runs out of products at the end of the planning horizon.This assumption renders stock-outs unavoidable, and thus restricts the applicability of the model.The introduced variant, referred to as the Cyclic Production routing Problem is aimed at producing practical solutions that can be repeated in cycles.Towards this direction.the ending inventory levels of nodes are considered equal to the starting ones.In addition, these inventory levels are treated as decision variables.As a result, the introduced CPRP model ensures that if no parameters are modified, the generated production and distribution plans are repeatable with cycle times equal to the length of the considered planning horizon more accurately reflecting real-world activities.Even though the cyclicity feature has been studied for various routing problems, to the best of our knowledge, this is the first time that it is incorporated and studied on the PRP model.The model is formulated as a mixed integer linear programming (MILP) problem.The formulation is strengthened with new and adapted valid inequalities of the PRP literature.A nonvehicle index flow formulation is proposed, thus the new CPRP model can be directly solved via a Branch-and-Cut algorithm (BnC) by any mathematical programming solver.
We compare the solutions obtained for the newly introduced cyclic and the basic PRP models for two wellknown data sets.To enable a fair comparison, the basic PRP is solved with starting inventories equal to the ones of the CPRP solutions.Due to the characteristics of the first data set, we also study a more applicable CPRP setting where the number of vehicles is restricted.Computational results indicate that the cyclic variant is significantly harder to solve: larger lower-upper bound gaps are observed and more computational time is required.When a restricted number of available vehicles is considered, the solution process becomes even harder, especially for larger instances.The total costs of the CPRP solution are from 25% up to 36%, and from 19% up to 22% higher than the PRP solutions, for the two tested data sets.The most important cost increase is seen for the variable production and inventory costs.In addition, the cyclicity requirement calls for the production and distribution of increased product quantities leading to significantly higher inventory levels.In terms of the vehicle fleet size, the cyclic variants require more routes compared to the classic PRP.This underlines that the basic PRP solutions underestimate the actual costs and cannot be directly adopted in real-life business settings, where the operational horizon is not finite.It is worth noting that the more applicable limited fleet configuration leads to higher inventory costs, but not considerably higher routing costs compared to the unrestricted vehicle model version.
Regarding future research directions, it would be interesting to investigate the effect of the length of the planning horizon.Also, a variant of the problem where each customer and route can have different cycles would be interesting to explore.Moreover, it is important to develop efficient algorithms capable of handling cyclicity, to enable solutions for large scale instances.In specific, such methodologies can either be local searchbased metaheuristic algorithms, or hybrid matheuristics incorporating mathematical programming components to handle the continuous quantities.A decomposition approach where the production and the visit decisions are taken in the first place, and routing is decided afterwards may also be effective against larger instances.
number of instances solved, O: number of instances solved optimally, I: number of instances that are infeasible, NS: number of instances not solved within two hours, Gap: average % gap of distance of the upper bound from the lower bound, T: average time in seconds.

Figure 3 .
Figure 3.Effect of cyclicity on the cost categories (data set A) (a) Set A1 (b) Set A2 (c) Set A3.

Table 2 .
Model comparison (data set A).

Table 3 .
Model comparison (data set B).

Table 4 .
Individual costs contribution to the total cost (%).

Table 5 .
Average inventory levels and total units stored (data set A).