An Inventory-Routing Problem with Pickups and Deliveries Arising in the Replenishment of Automated Teller Machines

T he purpose of this paper is to introduce, model, and solve a rich multiperiod inventory-routing problem with pickups and deliveries motivated by the replenishment of automated teller machines in the Netherlands. Commodities can be brought to and from the depot, as well as being exchanged among customers to efﬁciently manage their inventory shortages and surpluses. A single customer can both provide and receive commodities at different periods, since its demand changes dynamically throughout the planning horizon and can be either positive or negative. In the case study, new technology provides these machines with the additional functionality of receiving deposits and reissuing banknotes to subsequent customers. We ﬁrst formulate the problem as a very large-scale mixed-integer linear programming model. Given the size and complexity of the problem, we ﬁrst decompose it into several more manageable subproblems by means of a clustering procedure, and we further simplify the subproblems by ﬁxing some variables. The resulting subproblems are strengthened through the generation of valid inequalities and solved by branch and cut. We assess the performance of the proposed solution methodology through extensive computational experiments using real data. The results show that we are able to obtain good lower and upper bounds for this new and challenging practical problem.


Introduction
The purpose of this paper is to introduce, model, and solve an inventory-routing problem with pickups and deliveries (IRPPD) arising in the replenishment of automated teller machines (ATMs).Our study is motivated by the problem faced by a transporter responsible for such operations in the Netherlands.As in other cash-intensive economies, the Dutch credit institutions are gradually replacing regular ATMs by recirculation ATMs (RATMs), which are capable of accepting and dispensing banknotes, as well as checking their quality and authenticity.RATMs therefore provide customers the capability of both depositing and withdrawing cash.These machines provide tangible benefits to banks and customers, and are becoming the new standard in many markets.By the end of 2013, a total of 433,780 RATMs were in operation worldwide, corresponding to 17% of all ATMs (RBR 2014).The number of RATMs increased by 41% between 2012 and 2014, whereas the total number of ATMs grew by only 14% over the same period.RBR (2014) expects that the fast growth of the number and share of RATMs will continue to intensify.RATMs allow retailers to make deposits and enable customers to later withdraw the same cash.In this sense, RATMs are considerably more self-sufficient than the regular machines, which can only dispense cash.RATMs only require a visit to prevent the inventory level from reaching zero or from attaining the holding capacity of the machines when the supply and demand 1078 Transportation Science 50(3), pp.1077 -1091- , © 2016 INFORMS INFORMS are out of sync.When an RATM is empty, a customer can only use it to deposit cash, and when it is full only withdrawals are allowed.
In the problem considered here, the RATMs are replenished or partly emptied by using a fleet of rented armored trucks, which are very expensive.These trucks deliver cash from a depot to some machines, collect cash from some others to bring it back to the depot, or transfer cash between machines.The last operation reduces the routing cost and sometimes allows smaller or fewer trucks to be used.An important feature of the problem is the presence of inventory holding costs.Indeed, cash lying in a machine generates an implicit holding cost partly because it is insured and also because it incurs lost interest income.
To keep RATMs fully operational, that is when customers can use the RATM to both deposit and withdraw cash, one must solve an IRPPD.The IRPPD combines the features of two well-known classes of the vehicle routing problem: the inventory-routing problem (IRP) and the pickup and delivery problem (PDP), which we now briefly review.
The IRP belongs to the broader field of vendormanaged inventory systems in which a supplier coordinates the inventory management of a number of locations (Coelho, Cordeau, and Laporte 2014).In IRPs, the supplier must simultaneously decide when to visit its inventory locations, how much to deliver to each of them, and how to combine the deliveries into vehicle routes.There exist numerous variants of this problem (Andersson et al. 2010;Coelho, Cordeau, and Laporte 2014) and the related literature is rapidly expanding.The first exact algorithm for the singlevehicle IRP was based on branch and cut (Archetti et al. 2007).It could solve instances with up to 50 customers and three periods, and 30 customers and six periods.Archetti et al. (2012) later presented a hybrid tabu search matheuristic algorithm capable of dealing with larger instances and yielding solutions with very low optimality gaps.Coelho, Cordeau, and Laporte (2012) presented an adaptive large neighborhood search heuristic for the multivehicle IRP, and Coelho and Laporte (2013a) were the first to solve the multivehicle IRP exactly.Recently, multivehicle and multicommodity IRPs were also solved to optimality by Coelho and Laporte (2013b).All of these algorithms are based on branch and cut.For a recent survey of models and algorithms, see Coelho, Cordeau, and Laporte (2014).A related problem appears in maritime transportation, in which ships have to visit several ports to continuously deliver and pickup merchandise and commodities.Applications include the distribution of cement (Christiansen et al. 2011), chemical products (Dauzère-Pérès et al. 2007), and liquefied gases (Rakke et al. 2011), among others.For reviews of maritime transportation, see Christiansen, Fagerholt, and Ronen (2004) and Christiansen et al. (2013).
The PDP concerns the collection and distribution of one or several commodities to and from a set of locations.Berbeglia et al. (2007) classify PDPs and distinguish between three different problem structures: many-to-many (M-M), one-to-many-to-one (1-M-1), and one-to-one (1-1).The M-M structure means that each commodity may have multiple origins and multiple destinations, and that each location may be the origin or destination of multiple commodities.In 1-M-1 problems, some commodities are picked up at the depot and transported to some locations, whereas other commodities are picked up at these locations and transported to the depot.The 1-1 structure refers to a context in which each commodity has a single origin and a single destination, like in dial-a-ride problems (Cordeau and Laporte 2007).
Two important classes of PDPs with an M-M structure are the Swapping Problem (SP) and the One-Commodity Pickup and Delivery Traveling Salesman Problem (1-PDTSP).In the SP, introduced by Anily and Hassin (1992), each vertex of a graph provides a commodity and requests a commodity, possibly the same one.The problem is to design a least cost vehicle route to satisfy all requests.This problem is NP-hard on general graphs, but polynomial on some special structures (Anily, Gendreau, and Laporte 2011).Erdo gan, Laporte, and Cordeau (2010) have developed heuristics and a branch-and-cut algorithm for the multivehicle case.In the 1-PDTSP, each vertex either provides or requests a given amount of a single commodity, and a single vehicle route must be designed to transfer the right amounts of the commodity among vertices.The problem was introduced by Hernández-Pérez and Salazar-González (2003) and was solved by branch and cut (Hernández-Pérez and Salazar-González 2003, 2004a, 2007) and by heuristics (Hernández-Pérez and Salazar-González 2004b;Zhao et al. 2009).The 1-PDTSP arises in the rebalancing operations in shared bicycle systems (Benchimol et al. 2011;Chemla, Meunier, and Wolfler Calvo 2013;Contardo, Morency, and Rousseau 2012;Erdo gan et al. 2012;Erdo gan, Laporte, and Wolfler Calvo 2014;Raviv, Tzur, and Forma 2013).Some algorithms (Contardo, Morency, and Rousseau 2012;Dell'Amico et al. 2014;Raviv, Tzur, and Forma 2013;Shu et al. 2013) are capable of handling the multivehicle case.The problems encountered in RATM replenishment and bicycle repositioning are similar in the sense that they both consist of moving a commodity between locations to keep its level within given limits.The main differences lie in inventory holding costs, which are larger in the case of RATMs, and also in the planning horizon.Indeed, in the case of RATMs, planning is typically done over several days whereas in the case of bicycles, repositioning operations occur on a daily basis, sometimes only 1079 once during the night and often several times during the same day.
The PDP with a 1-M-1 structure deals with two types of commodities.Delivery commodities are transported from a single depot to multiple nodes, and pickup commodities are transported from multiple nodes back to the depot.A typical application is the recycling of products such as beer bottles, pallets, and containers (Berbeglia et al. 2007).A distinction is made between combined and single demands.Combined demands occur when locations may require both a pickup and a delivery, whereas single demand refers to problems where each location has either a pickup or delivery demand.The single demand problem was introduced by Mosheiov (1994).Heuristics were proposed by Gendreau, Laporte, and Vigo (1999) and by Subramanian and Battarra (2013), whereas Baldacci, Hadjiconstantinou, and Mingozzi (2003) and Hernández-Pérez and Salazar-González (2007) solved the problem exactly by branch and cut.The combined demand case was introduced by Min (1989).It was solved heuristically by Bianchessi and Righini (2007); Subramanian et al. (2010); Zachariadis, Tarantilis, and Kiranoudis (2010); Vidal et al. (2013), and exactly by branch-and-price algorithms by Angelelli and Mansini (2002); Dell 'Amico, Righini, and Salani (2006); Subramanian and Battarra (2013).We are not aware of any contributions on the multicommodity 1-M-1 PDP.
Our aim is to model and solve the IRPPD arising in the replenishment of RATMs.We view it as a novel logistics application arising in the banking service sector.As will be seen, the Dutch application that motivates this study gives rise to a very largescale mathematical model, which cannot be solved directly.To tackle it, we first decompose it into several subproblems through the application of a clustering procedure, and we further simplify the subproblems by fixing some variables.The resulting subproblems are then strengthened through the generation of valid inequalities and solved by branch and cut.The paper makes four main contributions.First, it introduces pickups and deliveries within an IRP context.Second, it combines two PDP structures: the 1-M-1 structure, which accounts for commodity movements from the depot to RATMs to the depot, and the M-M structure, which refers to commodity transfers among RATMs.Hence, our PDP could appropriately be designated as a 1-M-M-1 problem.Note that we tackle a rather general case of this problem, in that there are several vehicles and side constraints in addition to the standard IRP and PDP features.Third, to cope with realistic sized instances, we propose a practical decomposition methodology combining clustering, variable fixing, and branch and cut.Our fourth contribution is to apply our algorithm to data derived from a real-world case arising in the Netherlands.
The remainder of this paper is organized as follows.We formally describe the problem in mathematical terms in §2, where we provide a mixed-integer formulation.The algorithm is described in §3.This is followed by the results of an extensive computational study in §4 and by conclusions in §5.

Mathematical Programming Formulation
The IRPPD is defined on a directed graph = .The set of vertices is partitioned into , where is a set of depots and is the set of RATM locations.The set = i j i j ∈ i = j is the arc set.Each RATM incurs unit inventory holding costs i per period (i ∈ ), and has an inventory holding capacity C i .A handling cost is incurred at each depot, and is proportional to the quantity picked up and delivered at the depot.The length of the planning horizon is p, with discrete time periods t ∈ = 1 p .A set of rented armored vehicles k ∈ = 1 K , each with capacity Q k and average speed s k , is given.A renting cost k per period is incurred if vehicle k is used.Each vehicle is able to perform one route per period, from one depot to a subset of RATMs, each requiring r units of time to be served, and back to the same depot.The shift of each vehicle is limited to S minutes, after which monetary units per minute of overtime are incurred.A travel time c ij in minutes is associated with arc i j ∈ .We assume each depot has sufficient inventory and capacity to perform all pickups and deliveries during the planning horizon.The inventories are not allowed to exceed the holding capacity nor are they allowed to become negative.At the beginning of the planning horizon, the decision maker knows the current inventory level I 0 i of the RATMs and receives information on the net demand d t i of each RATM i for each period t.Negative demands mean that the RATM provides a commodity, whereas positive demands mean that the RATM receives some quantity of the commodity.We assume that the quantities q t i received by RATM i in period t can be used to satisfy its net demand in that period.The quantities picked up at the depots and at the RATMs may be delivered to any RATM to satisfy their demands.The objective of the problem is to minimize the total cost while satisfying the net demand for each RATM in each period.
The variables used in the formulation are as follows.The multidepot IRPPD is then formulated as follows: subject to The objective function (1) minimizes the total cost of inventory holding, inventory handling at the depot, and vehicle renting.Constraints (2) state the inventory conservation condition over successive periods: they define the inventory in period t as the inventory held in period t − 1, plus the quantity delivered, minus the quantity picked up, plus the net demand of the location.Constraints (3) define the bounds on the inventory held by each RATM throughout all periods.Constraints (4)-( 9) guarantee that proper vehicle routes are created.More specifically, constraints (4) are flow conservation constraints, constraints (5) relate the incoming flow of a vertex to the fact that it is visited or not.Constraints (6) state that a vehicle can leave at most one depot in any period.Constraints (7) mean that a vertex cannot be visited if no vehicle is assigned to it.Similarly, constraints (8) relate the use of a vehicle during a given period to the flow variables associated with that vehicle and that period.Constraints (9) prevent the formation of subtours over a set of RATMs.These are the subtour elimination constraints proposed by Gendreau, Laporte, and Semet (1997) for the covering tour problem.Constraints (10) link the delivery binary variables w kt i to the visiting binary variables y kt i .They allow a delivery decision to be made only if a visit is performed to the RATM.Constraints (11) are similar to (10) and apply to the pickup decisions.Constraints ( 12) and ( 13) allow a quantity to be delivered or picked up at an RATM only if the corresponding binary decision variable is equal to 1. Constraints ( 14) ensure that a visit to an RATM is made either for a pickup or for a delivery operation.Constraints (15) guarantee that the shift duration is respected and that overtime is properly accounted for.Constraints (16) ensure that the load within the vehicle is consistent along its route.Constraints (17) set bounds on the load inside the vehicle after serving RATM i. Constraints (18) set the load of the vehicle when leaving the depot, whereas constraints (19) compute its load when returning to the depot at the end of the route.Finally, constraints (20)-( 23) define nonnegativity and binary conditions on the variables.

Solution Algorithm
In the RATM application that motivates this study, 6,377 machines must be replenished by 32 vehicles, 1081 each based at a different depot, over a six-day planning horizon.This means that the (MD-IRPPD) model just described will contain more than 7.8 billion binary variables and about the same number of nonlinear constraints, in addition to O 2 subtour elimination constraints.Solving such a model or even its continuous relaxation is beyond the reach of any available mathematical programming solver.We have therefore opted for a practical three-step solution strategy.We first apply a clustering procedure to decompose the initial problem into 32 subproblems of about 200 RATMs each, one for each depot (each having a single vehicle).This reduces the number of variables of each subproblem to a more manageable size, around 250,000, but still too many to reach optimality.For comparison purposes, the largest IRP instance that can be solved exactly, with a single vehicle and a single depot, contains 200 customers (Coelho and Laporte 2014), but the IRP does not contain many of the features presented in this paper.Likewise, the largest M-M PDP instance solved to optimality contains 200 vertices and a single vehicle (Hernández-Pérez and Salazar-González 2007).The IRPPD combines features of these two problems and is significantly more complicated than either of these.For these reasons, it is unrealistic to solve this problem exactly for the large sizes we are considering.We therefore simplify the subproblems further by fixing some variables.This yields instances that can be solved exactly or to near optimality by branch and cut.These three phases are described in § §3.1-3.3,respectively.Algorithm 1 provides a sketch of our solution procedure.
Algorithm 1 (Overall solution algorithm) 1: Clustering procedure: solve a transportation problem to obtain as many single-depot instances as the number of depots (see §3.1) 2: for each single-depot instance do 3: Apply the variable fixing procedure (see §3.2) to reduce the instance size 4: Solve the IRPPD model by branch and cut with some fixed variables and valid inequalities (see §3.3) 5: end for.

Clustering Procedure
Partitioning the customer base into clusters not only makes sense from a computational point of view but also from an empirical standpoint because this is how the problem is solved in practice.We solve a transportation problem to assign RATMs to depots while controlling the size of the clusters.The assignment cost in each cluster is the sum of the radial average travel times, interpreted here as distances.Eilon, Watson-Gandy, and Christofides (1971) show that radial distances are often a good proxy for routing costs.Let ij be a binary variable equal to 1 if and only if RATM i is assigned to depot j.The transportation problem is then minimize The objective function ( 24) minimizes the total sum of the distances between assigned RATMs and depots.Constraints (25) ensure that each RATM is assigned to exactly one depot, and constraints (26) ensure that each depot is assigned to at most A RATMs, with A chosen so as to balance the number of RATMs per depot.Constraints ( 27) define the domain of the variables.
After the clustering procedure has been applied, the problem is simplified and the (MD-IRPPD) formulation can be reduced to one without the depot index, i.e., = 0 .We provide this formulation for the sake of completeness and because we will later use it to generate valid inequalities.Reusing the notation employed in §2, we now state the single-depot formulation subject to The description of the variables, objective function, and constraints is similar to that of the multidepot model.(2009).When it is known in advance whether an RATM requires a pickup or a delivery, they can be lifted (Gribkovskaia et al. 2007).However, this is not the case here.The third comment refers to the fact that the formulation is nonlinear because of constraints ( 36), (37), (40), and (43).However, these can be linearized as follows.Constraints ( 36) and ( 37) can be rewritten in a slightly weaker form as On their own, these constraints would allow quantities to be picked up or delivered to violate the inventory bounds, but together with constraints ( 29) and (30), they are feasible linear representations of constraints ( 36) and (37).Constraints (40) can be linearized as Finally, constraints (43) can be rewritten in a linear form as The IRPPD formulation can be strengthened through the generation of the following valid inequalities: Constraints ( 53) are referred to as logical inequalities.They tighten the relation between routing and visiting variables.Constraints (54) include the supplier in the route of vehicle k if any RATM is visited by that vehicle in that period.

Variable Fixing Procedure
To solve realistic instances, we have designed a variable fixing procedure, which is applied prior to executing the branch-and-cut algorithm, in the spirit of granular search (Johnson and McGeoch 1997;Toth and Vigo 2003).This heuristic determines for each RATM i ∈ and each period t ∈ whether it must be visited, must not be visited, or must perhaps be visited.Specifically, for each period t we distinguish five subsets: (1) a must pickup set t , (2) a must deliver set t , (3) a perhaps pickup set t , (4) a perhaps deliver set t , and (5) a not visit set t .These five sets define a partition of .When RATMs are added to one of the subsets t , t , or t , the size of the formulation is considerably reduced and the branch-and-cut algorithm benefits from this reduction, because some of the decisions related to pickups, deliveries, and visits for some RATMs are already made.When RATMs are added to one of the other two subsets t or t , the algorithm still determines for each period t whether a visit is required, but the benefit of the variable fixing procedure is that it limits the choice to either a pickup or a delivery.It is implemented as follows: Also, it follows directly that several other constraints fixing binary variables can be derived from ( 55)-( 59).These are used to further reduce the size of the problem by effectively eliminating several variables from the problem k∈ y kt i = 1 i ∈ t t ∈ (60) The variable fixing procedure uses the parameters defined below, which allow us to test various settings.The first parameter relates to due periods and tipping periods.A due period is a period in which a visit is required to prevent the inventory from exceeding the holding capacity or from a stock out.A tipping period is a period in which a pickup is required to remain cost efficient, that is when the holding cost is disproportionately high related to the cost of visiting the RATM.The parameters are as follows: • m: number of periods preceding a due period or a tipping period; • f : number of RATMs that are closest to RATM i; • g: minimum inventory level of RATM i, in a percentage of its holding capacity C i ; • b: expected number of RATM visits per vehicle k in a period t.
In what follows, let Īt i = I 0 i + t t=1 d t i be the cumulative inventory of RATM i in period t .This is useful to identify up to which period t the RATM will be capable of respecting its inventory constraints without replenishment.The pseudocode of the variable fixing procedure is presented in Algorithm 2. The variable fixing procedure starts by adding each RATMs i in each period t to the not visit subset t (lines 2-4).When t = 1 or m = 0, lines 5-10 ensure that RATMs are added to the must visit subsets t and t when the holding capacity C i is exceeded or when the inventory Īt i becomes negative.If either condition is satisfied, the algorithm does not have the flexibility to choose between preceding or succeeding periods to visit the RATM.Therefore, the RATM must be visited.
If the condition in line 5 is not met, we check three more conditions in lines 13-18 to add RATMs i to the perhaps visit subsets t and t in period t as well as to the subsets in m preceding periods t.Line 12 ensures that RATMs i are also added to the subsets in m number of preceding periods t.We introduce lines 13 and 14 to ensure RATM i is visited for a pickup when a tipping period occurs, i.e., when the holding cost exceeds the cost of visiting the RATM.More precisely, the following condition is verified: the inventory Īt i of RATM i in period t multiplied with the unit inventory holding cost i exceeds the vehicle rental cost k divided by the expected number of visits per route b and divided by the expected number of days elapsed since the last visit.To this end, we calculate the expected number of days elapsed since the last visit by dividing the inventory level Īt i by the average demand t∈ d t i /t.If the latter condition is met, then the RATM i is added to the perhaps pickup set t .RATMs i are also added to the t subset in periods t when the inventory exceeds capacity Īt i > C i (lines 15 and 16).In the case of stock outs, the RATMs are added to the perhaps deliver subset t , which is indicated in lines 17 and 18.
Most likely, some RATMs i will have been taken out of set t and added to the other subsets t t t , and t when arriving at line 18.The final part in lines 22-28 ensures that even more RATMs j, which are located closest to RATMs i belonging to t t t , or t , are added to the perhaps visit subsets.The number of RATMs j added per RATM i is determined 1084 Transportation Science 50(3), pp.1077-1091, © 2016 INFORMS by parameter f .It should be noted that RATMs j are only added when they meet an inventory level Īt j , which is at least as high as a given percentage g of the holding capacity C i for a pickup, or lower than a given percentage 1 − g for a delivery.The rationale of this final part is to stimulate the exchange of inventories between RATMs i and j to eventually reduce the amount to be picked up from, and delivered to the depot.

Branch-and-Cut Algorithm
We have implemented a branch-and-cut algorithm capable of solving the IRPPD model with the linearized constraints, the valid inequalities, and the fixed variables.All variables of the formulation are explicitly handled by the algorithm, but we cannot generate all subtour elimination constraints (33) a priori.These will be dynamically generated as cuts as they are found to be violated.The formulation is then solved by branch and cut as follows.At a generic node of the search tree, a linear program with relaxed integrality constraints is solved, a search for violated constraints is performed, and violated valid inequalities are added to the current program, which is then reoptimized.This process is reiterated until a feasible or dominated solution has been reached, or until no more cuts can be added.At this point, branching on a fractional variable occurs.Details regarding the implementation and improvements to this algorithm are provided in §4.2.

Computational Experiments
We will present the instances generator in §4.1, some implementation details in §4.2, the analysis of our extensive computational experiments in §4.3, and the estimated benefits of our study in §4.4.The instance set as well as detailed computational results are available at http://www.leandro-coelho.com/instances.

Instances Description
The data used in our experiments stem from a realworld case in the Netherlands.A total of 6,377 cash dispensing self-service devices, both regular and recirculation ATMs, were installed in the Netherlands in 2013 (ABN ARMO Bank 2013;ING Bank 2013;Rabobank 2013).These are replenished from 32 cash centers (depots).By setting A equal to 200 in constraints (26), we assigned the 6,377 RATMs fairly evenly to the depots, each of which serves a small region of the Netherlands (Figure 1).For our experiments, we assume that all self-service devices in the Netherlands are RATMs.This is not unrealistic given the recent strong increase in the share of RATMs in the Netherlands (ECB 2013).Also, in several other countries, such as Japan, RATMs already dominate the cash self-service device market (RBR 2014).To conduct our numerical experiments, we used real data for the locations of RATM devices, true distances and speeds, and real RATM demands.Although the Netherlands is rather densely populated, not all parts of it are equally well served by RATMs.Their coverage varies considerably over the country.Our solution methodology is rather robust and can easily deal with this diversity.Figure 2 depicts the position of the RATMs and the cash center for the Amsterdam area.There still exist some ARMs that are not RATMs.For these, we generated RATM demands.The depot locations are confidential for security reasons, so we used approximate locations.
To mimic real demands for devices that are not yet RATMs, we have carefully studied real demands from readily installed RATMs.This study demonstrates that the number of transactions per day is well estimated by a Poisson probability distribution; the number of deposits and the number of withdrawals follow completely different patterns in terms of quantity, value, and weekly seasonality pattern; a huge diversity exists among the demand intensity of RATMs, and the average value per deposit is E1,000 and the average value per withdrawal is E133.
Based on these findings, we have generated deposit and withdrawal distributions separately and added these sequentially to construct a net demand per RATM per day.For the deposits and withdrawals of each RATM, we uniformly assigned two demand intensity levels of one to five to mimic demand intensity for both deposits and withdrawals (Table 1).
By randomly drawing from a Poisson distribution with the assigned number of transactions per day and multiplying the resulting number of transactions by the average transaction amount, i.e., E1,000 or E133, we obtained daily deposit and withdrawal amounts.Finally, to cope with week seasonality, we modified these amounts with the factors described in Table 2.
The parameters of our instance are then as follows: • i = E0 08 per E1,000 inventory per period t (based on a 3% annual interest rate).
• C i = E260,000 per RATM i.This is an estimation based on an RATM with four cassettes, each with a capacity of 2,000 notes.
• p = six periods.Periods coincide with days, which together define a workweek from Monday to Saturday.In practice, order lead times are not more than one or two days and so a planning horizon of six days is sufficient.
• = one vehicle.• Q k = E7,800,000 per vehicle k.This is an estimation based on 30 full replenishments of E260,000 in a single period.• k = E2,000 per vehicle k per period used.
• S = eight hours.This is the regular daily work time in the Netherlands.We assume the vehicle is loaded prior to performing the route, so the driver has at most eight hours to perform all pickups and deliveries.
• = E8 00 per minute, which is approximately twice as expensive as the regular vehicle renting cost.
• I 0 i ∈ 0 C i .Each RATM i is assigned a random initial inventory in euros.
• d t i = ± − 45,000 ±45,000 per period.Random values are drawn from a Poisson distribution for both withdrawals and deposits, which are thereafter combined into a net demand.To simulate real demands at different locations, we have ensured that RATMs either have a net-positive, net-negative, or balanced demand without losing stochasticity.We have then generated different demands at distinct periods: in periods t ∈ 1 2 the demand tends to be net-positive, because in the Netherlands more cash is deposited on Mondays and Tuesdays.The demand in periods t ∈ 5 6 tends to be net-negative, since more cash is withdrawn on Fridays and Saturdays.
• c ij = the arc set is constructed by first using real travel distances between the RATMs and the depot.A routable network data set for the Netherlands was then constructed using OpenStreetMap data (Geofabrik GmbH, OpenStreetMap Contributors 2013) and average driving speeds were used on the various road types to approximate true speeds.

Implementation Features
The algorithm just described was coded in C++ using IBM Concert Technology and solved with CPLEX 12.5.1 running on a single thread.All computations were executed on a grid of Intel Xeon processors running at 2.66 GHz with up to 48 GB RAM installed per node, with the Scientific Linux 6.1 operating system.A time limit of six hours was imposed on the execution of each of the 32 instances.
The clustering procedure was run only once.It is conceivable that overall improvements could have been obtained by performing 1-opt or 2-opt moves between clusters, but there exist O 2 such potential moves and evaluating each of them would have entailed solving a large-scale integer linear program, which was impractical.
Several experiments were performed with various variable fixing parameters.These settings were gathered from cash supply chain parties in the Netherlands, and estimated when these could not be made publicly available for security reasons.The parameters for the clustering heuristic were set as follows: m ∈ 0 1 2 ; f ∈ 0 1 2 ; g ∈ 30 50 70 ; b ∈ 15 .
For the instances considered in this paper, which contain up to 200 vertices, we cannot obtain optimal solutions with our branch-and-cut algorithm.To cope with this situation, we have decided to add two new layers to the branch-and-cut algorithm.In the first one, we verify at every node having a fractional solution whether inequalities (53) are violated and we then add them to the linear program.This helps improve the lower bound of the instance.The second new layer is added to handle constraints (50) since we have observed that these constraints are not tight with respect to the vehicle capacity and are generally satisfied.For this reason, we have devised a procedure to verify whether they are violated only at an integer solution.If the integer solution is found to violate constraints (50), these are added and the node is then reoptimized.Otherwise, the solution is feasible for the IRPPD.These two enhancements have significantly decreased the computational time.
Moreover, we have also observed that the separation algorithm of constraints (9) can be improved by adding a new set of variables representing the route of the vehicle in an undirected graph.By doing this, the separation algorithm runs on a graph half the size of the original one, looking for connected components and deriving maximum cuts over a much smaller network.Moreover, each subtour elimination cut derived for this new variable is equivalent to two cuts expressed in the original variables, one in each direction.
These two procedures employed to generate new cuts dynamically typically yield around 20,000 cuts only at the root node of the search tree.This number, although sizeable, is only a small fraction of the total number of cuts that could potentially be generated.By optimizing a problem with fewer constraints, we gain in speed since far fewer simplex iterations are needed.We note that a typical IRPPD instance still contains around 250,000 binary variables and 150,000 constraints after the two procedures just described have been applied, besides O 2 subtour elimination constraints.
We have also observed that at the beginning of the optimization process, the problem is degenerate, i.e., several pivoting operations do not improve the value of the objective function.We have therefore taken advantage of the fact that during the optimization process, one can change the routing decisions while remaining feasible.Since routing variables x do not appear in the objective function, this solution has the same objective function value.Note that as long as the total route length is less than the shift duration, the solution remains feasible.We have also tested adding a small coefficient to the objective function to further minimize route length, and thus avoid degeneracy, but this option did not yield any significant advantage.
Finally, to obtain a clear view of the performance of our algorithms, and also to speed up the solution of each node, we have turned off the CPLEX cut generation.

Analysis of the Computational Experiments
We start our analysis by presenting the results obtained by the algorithm presented in § §3 and 4.2 on the set of 32 instances described in §4.1.Not all conceivable instances of the MD-IRPPD are feasible, but feasibility was always achieved with our data.The optimality gaps observed before variable fixing are rather large, with an average of 51%, a maximum of 62%, and a minimum of 45%.Even when 24 hours of computing time were allotted, these figures did not improve significantly: the average lower bound increased by 2% from 16,002 to 16,355, and the average gap went down from 51.31% to 50.23%.In what follows, we use the information obtained within six hours of computing time to allow for a direct comparison with the performance of the variable fixing procedure.To benchmark the quality of these solutions, we can compare to the instances of the IRP without pickup and delivery, which has recently been solved to optimality in Coelho and Laporte (2014) for comparable sizes.The problem at hand is much more complicated and clearly requires an additional effort to obtain good solutions.
We have then limited the flexibility of the algorithm by disallowing visits to RATMs that do not require a pickup or delivery to remain operational.This is achieved by setting to zero both parameters m and f of the clustering procedure.Obviously, all solutions 1090 Transportation Science 50(3), pp.1077-1091, © 2016 INFORMS being transported, the more expensive the insurance becomes.Our algorithm achieves a 54% reduction in the in-transit inventory, when compared with the traditional 1-M-1 structure.Not only will the insurance premium decrease substantially, but a less expensive truck (due to reduced armor) could be deployed, CO 2 emissions can be reduced, and losses due to thefts should significantly go down.

Conclusions
We have introduced, modeled, and solved an inventoryrouting problem with pickups and deliveries.We have been successful in solving a difficult application of the problem arising in the optimization of distribution and inventory management of cash in recirculation ATMs in the Netherlands.The problem was first modeled as a very large-scale nonlinear mixed-integer program.After applying a clustering procedure and linearizing some of the constraints, including some valid inequalities and fixing the values of several variables, we have solved the problem by branch and cut.Different settings of the variable fixing procedure were tested, and we have shown which ones are able to provide a good trade-off in terms of simplification of the problem and upper bound values.After variable fixing, we were able to solve exactly or to near optimality 32 instances involving up to 200 vertices.This size is similar to that of the largest IRP or M-M PDP instances that have been solved in the past.Our problem is obviously more difficult, because it combines these two features.Since recirculation ATMs are expected to become a market standard in the near future, our paper anticipates this development by providing cash supply chain parties the means to immediately improve the replenishment operations and yield significant cost savings. 1085 Figure 1 (Color online) Map of the Netherlands Depicting 32 Subregions

Figure 2 (
Figure 2 (Color online) Map of the Amsterdam Subregion Five families of binary variables are used: directed routing variables x kt ij are equal to 1 if and only if arc i j is used on the route of vehicle k in period t; visiting variables y kt Transportation Science 50(3), pp.1077-1091, © 2016 INFORMS period t; and pickup variables z kt i are 1 if and only if a pickup is performed at RATM i by vehicle k in period t.Integer variables I t i represent the inventory level at RATM i ∈ at the end of period t ∈ .Variables q kt i are integers representing the product quantity delivered to RATM i using vehicle k in period t, and variables p kt i are integers representing the product quantity picked up from RATM i using vehicle k in period t.Two sets of variables represent the amount of inventory carried by vehicle k in period t into and out of the depot: H t k and J t k , respectively.The load of vehicle k after serving RATM i in period t is represented by variables u kt i , and E t k represents the overtime in number of extra minutes worked by vehicle k in period t beyond the regulatory shift duration S.
i are equal to 1 if and only if RATM i is visited by vehicle k in period t; vehicle usage variables v kt are equal to 1 if and only if vehicle k is used in period t; delivery variables w kt i are equal to 1 if and only if a delivery is made to RATM i by vehicle k in 1080 Miller-Tucker-Zemlin (1960) with respect to the IRPPD formulation.The first regards constraints (40), which resemble theMiller-Tucker-Zemlin (1960)subtour elimination constraints.These constraints do not eliminate subtours in this context, because the load within the vehicle is not monotonically increasing or decreasing.We therefore imposeDantzig-Fulkerson- Gribkovskaia et al. (2007)ination constraints (33) whose number is O 2 .Second, the load consistency constraints (40) have been used in a number of other pickup and delivery problems; see, e.g.,Desaulniers et al. (2002);Gribkovskaia et al. (2007);Hoff et al.

Table 1
Different Levels of Deposits and Withdrawals per Day