Synchronizing vans and cargo bikes in a city distribution network

One of the significant side-effects of growing urbanization is the constantly increasing amount of freight transportation in cities. This is mainly performed by conventional vans and trucks and causes a variety of problems such as road congestion, noise nuisance and pollution. Yet delivering goods to residents is a necessity. Sustainable concepts of city distribution networks are one way of mitigating difficulties of freight services. In this paper we develop a two-echelon city distribution scheme with temporal and spatial synchronization between cargo bikes and vans. The resulting heuristic is based on a greedy randomized adaptive search procedure with path relinking. In our computational experiments we use artificial data as well as real-world data of the city of Vienna. Furthermore we compare three distribution policies. The results show the costs caused by temporal synchronization and can give companies decision-support in planning a sustainable city distribution concept.


Introduction
The World Health Organization states in its Annual Report 2013 that two of the major challenges in the upcoming years are the aging of the population and growing urbanization, given that 70 % of the world population is forecast to live in cities by 2050 (World Health Organization-Centre for Health Development 2013). Urbanization can have a great impact on population health and the environment and is therefore a keyword high on the agenda of authorities and the public.
One important aspect of urbanization is the increasing traffic volume in cities, and with it the potential rise in traffic congestion, pollution and noise nuisance. Ecommerce and home delivery services are favored distribution channels, required among others by the elderly population or persons with complex daily schedules who cannot easily make their retail purchases on-site.
Thus, accessibility and mobility for citizens and freight need to be improved in urban areas. A means of reducing inefficient polluting traffic within cities is to switch from conventional vehicles to other transport modes when entering urban zones. This of course requires coordination between different stakeholders in order to achieve an optimal and sustainable utilization of resources. Any re-organization of city logistics must be managed in such a way that traffic is reduced and the modal split shifts towards environmentally friendly transportation modes.
Especially in urban areas, there is a huge potential for consolidation and coordination of distribution flows, since typically small amounts of freight need to be delivered to spatially dispersed locations. The aim of the current paper is to look at how to efficiently organize the distribution of goods in cities by consolidating the transport requirements of different stakeholders and using environmentally friendly transport modes in inner-city areas. Reduced traffic volume, efficiently utilized transport resources and the use of green transport modes counteract negative consequences for the environment and the population and enable more livable urban areas in the future.
To identify the main problems for companies concerning city distribution of goods in Vienna we conducted expert interviews with logistics managers from two different sectors: pharmacy wholesale and distributors of vegetable boxes (veg boxes), who distribute not only locally-grown and organic vegetables but also fruits and sometimes even meat and dairy products from the grower or a small co-operative directly to the end consumer (Brown et al. 2009). We have chosen pharmacy wholesale because this sector requires a relatively high frequency of deliveries during the day. Each of the 316 pharmacies (status 2013) in Vienna is supplied at least twice a day (Österreichische Apothekerkammer 2014). The interviewed distributors of veg boxes, on the other hand, are looking for new and more sustainable distribution concepts, although they serve their customers-mainly end consumers-only once a week.
Nevertheless, there are a lot of similarities between those two sectors as our expert interviews have shown: Goods from a variety of suppliers are consolidated at a depot at the outskirts of the city and then delivered by small vans to the customers. Customers are often located at poorly accessible parts of the city such as in pedestrian zones or zones with access restrictions. The main difficulties for the companies result from narrow streets and missing parking space, especially in the city center, which often leads to congested streets. For the companies these circumstances result in vehicles that are not fully loaded, and therefore to inefficient operations. Thus, costs and emissions could be reduced by switching to a more efficient distribution policy.
In this paper we investigate three possible distribution policies to help solve these issues. The first policy reflects the current delivery scheme used by the pharmacies in which customers are supplied directly from the main depot. This corresponds to a classical vehicle routing problem (VRP), in which inner-city areas are supplied by trucks. This can involve difficulties since a lot of cities have already implemented access restrictions for conventional vehicles. In many cities, the city center is a place of historic interest with narrow, congested streets in which additional constraints on delivery apply such as restricted space, accessibility constraints or time windows for loading. Moreover, it is the heart of the city where tourists and inhabitants want to enjoy their time without having to cope with the nuisances of freight traffic.
To tackle these problems we look at a second distribution policy in which customers in the inner-city are supplied by smaller, eco-friendly vehicles. We investigate a twoechelon distribution scheme where vans perform the delivery outside the city center, on the so-called first echelon, and supply the eco-friendly vehicles, that operate on the so-called second echelon. The transfer of load between first and second echelon vehicles is performed on satellite facilities located outside the inner-city area. The downside of this policy is that space is rare in urban areas and hence possible satellite locations are difficult to acquire.
Finally, to overcome this drawback we investigate a third policy that uses mobile satellites without storage facilities, which, for example, could be a parking space. With this policy temporal synchronization of the first and second echelon vehicles at those satellites is necessary.
The contributions of this paper are the following. First, we develop a solution algorithm based on a greedy randomized adaptive search procedure (GRASP) with PR for a two-echelon, multi-trip vehicle routing problem with temporal synchronization. Our algorithm is fast and efficient, and it can easily be used by urban planners. Furthermore, we look at a case study of the city of Vienna and investigate the three different distribution policies explained above.
The paper is organized as follows. Section 2 provides a literature review. Section 3 describes the problem that we study in detail. Section 4 presents the solution method that we develop. Section 5 shows our computational results, while Sect. 6 concludes the paper.

Literature review
We consider a two-echelon, multi-trip vehicle routing problem with synchronization. The classic two-echelon vehicle routing problem (2eVRP) was first solved by Perboli et al. (2011). Since then, exact methods (Jepsen et al. 2013;Baldacci et al. 2013;Santos et al. 2014) and heuristics (Hemmelmayr et al. 2012;Crainic et al. 2011b) have been proposed. In the 2eVRP studied in these papers, satellites are assumed to have (little) storage space so that temporal synchronization is not necessary. Moreover, multiple trips are not included, and the second echelon vehicles are based at one satellite and return there at the end of the day.
Another aspect of the problem described in this paper is vehicle synchronization, i.e. vehicles of both echelons have to meet at a certain place (spatial synchronization) at a certain point in time (temporal synchronization). A survey of synchronization in vehicle routing is given in Drexl (2012), and a classification of different types of synchronization is provided.
For a detailed review about vehicle routing problems in city logistics we refer the reader to the survey of Cattaruzza et al. (2015). Furthermore, a review of two-echelon routing problems can be found in Cuda et al. (2015).
In the following, we will focus on those papers that consider problem aspects which are similar to those dealt with in this paper. Crainic et al. (2009) study a time-dependent, two-echelon, multi-trip routing problem with synchronization. However, no solution algorithm is proposed, and only promising directions are indicated.
In Crainic et al. (2011a) a modeling framework for tactical planning in city logistics with uncertain demand is developed. They propose a two-stage stochastic programming formulation with recourse. In the first stage the routing of customer demand to the synchronization points is determined. The recourse strategy of the second stage has to solve a synchronized, scheduled, multi-depot, multiple-tour, heterogeneous VRPTW. Nguyen et al. (2013) propose a tabu search for the time-dependent multi-zone multitrip vehicle routing problem with time windows. This problem occurs on the second echelon of two-echelon city distribution schemes. The arrival of first-echelon vehicles at synchronization points determines the availability periods in which the second-tier vehicles must arrive at satellites, load and then service customers. Grangier et al. (2014) study a two-echelon multiple-trip vehicle routing problem with satellite synchronization and time windows. They propose an adaptive large neighborhood search heuristic to solve the problem. Computational experiments are performed on adapted Solomon instances for the VRP. Though their problem is related to ours, they do not include visits to regular customers on the first echelon and they consider time windows at the customer nodes.

Problem description
The problem we seek to solve in this paper is a two-echelon routing problem with synchronization, in which the inner-city delivery on the second echelon is performed by cargo bikes. The customers are divided into two groups. There are those located within the city center called bike-customers, and those located outside this area. These latter customers are supplied by vans and hence referred to as van-customers. Moreover, the vans also supply goods to the cargo bikes via so-called satellites. The satellites are transshipment points where vans and cargo bikes can meet. After loading, the cargo bikes perform their delivery and when they have to reload, they move again to a satellite. The city center is defined as the area within the satellites. As we want to avoid vans completely in the city center we also penalize van trips crossing this area by imposing penalty costs for such trips.
Goods from different suppliers arrive at the van depot located on the outskirts of the city. From there, vans perform the delivery to van-customers and to satellites. The satellites are a set of predefined meeting points which do not have storage facilities.
So vans and cargo bikes must meet in a synchronized way at the same time at the same physical satellite, while waiting times for cargo-bikes and vans are minimized. Similar to vans, cargo bikes start and end their tours at the bike depot located in the city center.
The formulation of our VRP is based on the ones by Grangier et al. (2014) and by Crainic et al. (2012). The problem is defined on a graph G = (V, A). The set of nodes V consists of the set of two depot nodes V d (one as starting point for bikes and one as starting point for vans), the set of n customers V c , and the set of m satellites V s .
Based on the topography of the problem the set of customers V c is split into two parts: All customers located in a circle C I ⊆ P V s (the polygon spanned by the satellites s ∈ V s ) are served by cargo bikes. We denote this subset of bike-customers by V b c and its complement V c \V b c by V v c (the van-customers). Let n 1 :=| V b c | and n 2 :=| V v c |. Furthermore, we can denote the bike depot by v b d , and the van depot by v v d . As each physical satellite s ∈ V s can be used as supply point for each bike-customer d ∈ V b c a set of cloned satellitess ∈Ṽ s is created where each physical satellite s is duplicated n 1 times (once for each bike-customer). So, the extended set of vertices can be denoted by The set of arcs A consists of all feasible arcs which is described in detail below. The first level (served by vans) is defined on has its specific loading or service time λ i ≥ 0 and each vertex i ∈ V c has its specific demand d i > 0.
Furthermore, we have two fleets of vehicles: a fleet of vans (F v ) and a fleet of cargo bikes (F b ) in non-limited number with F = F v ∪ F b . Each vehicle type κ = {v, b} has its specific speed s v /s b to derive the travel times , if now separate travel time matrix is provided, and its capacity Q v /Q b . For the calculation of costs each vehicle type has vehicle cost c v D /c b D per unit of distance traveled, driver cost c v T /c b T per unit of time used and fixed cost c v F /c b F per vehicle used. To penalize vans crossing the inner area C I ⊆ P V s we use a penalty cost matrix with penalty costs p i j for each arc (i, j) ∈ A v . Please note that p i j is 0, if an arc does not cross C I . So c k i j represents the costs of vehicle k for traversing arc (i, j). Furthermore, we define the following decision variables: The binary variable 1 if vehicle k travels from node i to node j 0 otherwise and the continuous variables -w b i ≥ 0 specifying the occurring waiting time for a bike at node i w v i ≥ 0 specifying the occurring waiting time for a van at node i t k i ≥ 0 specifying the arrival time of vehicle k at node i u b i ≥ 0 specifying the load of a bike after serving node i u v i ≥ 0 specifying the load of a van after serving node i Furthermore, we use the constant M u which represents the total demand of all customers.
The maximum route duration is given by t max and the maximum permitted waiting time at any satellites ∈Ṽ s for each vehicle is defined by w max .
The objective function of our model can then be formulated as: Subject to: The objective function (1) minimizes the total transportation costs, including variable costs based on distance and time traveled as well as fixed cost for each vehicle in use. Constraints (2) and (3) ensure that each van leaving the van depot also returns to this depot, while constraints (4) and (5) enforce same for bikes. Constraints (6) and (7) guarantee that each customer is visited exactly once by a vehicle. Constraints (8), (9), (10) and (11) ensure that each cloned satellite can be used at most once by a bike or a van and if it is used as supply point by a bike it has to be visited by a van too. Constraint (12) ensures for each bike that the load after visiting a node is equal to the load before minus the demand of the respective node. Constraint (13) ensures that each bike is empty when reaching a satellite. Constraint (14) ensures that the load of a bike at each visited node does not exceed the bike capacity. Constraint (15) ensures for each van that the load after visiting a node is equal to the load before minus the demand of the respective node. Constraint (16) ensures that the load of a bike at a satellite is added as demand for the van which visits this satellite. Constraint (17) ensures that the load of a van at each visited node does not exceed the van capacity. Constraints (18) and (19) restrict the maximum tour duration of a bike and a van, while constraint (20) enforces the start time of vans at the depot to be zero. Constraints (21) and (22) ensure the scheduling of nodes for bikes and vans. Constraint (23) determines the waiting time of a bike at a node as the difference between the arrival time of the van at this node minus the arrival time of the bike at this node, and if the difference is negative the waiting time of a bike at this node is set to zero. Constraint (24) ensures the same issue for the waiting time of a van at a node. Constraints (25) and (26) ensure the synchronization of bikes and vans at cloned satellites within a maximum permitted time span. Eventually, constraint (27) imposes binary values on the decision variable, while constraint (28) ensures non-negative values for all other variables.
The used notation is given in Table 1. Figure 1 provides a simplified representation of our routing problem which shows the van depot (large rectangle), the bike depot (small rectangle), a number of potential satellites (triangles) and a number of customers (dots). The customers within the inner circle C I of the satellites are served by cargo bikes (dashed lines) and all other customers are served by vans (solid lines). At a satellite a cargo bike has to meet with a van to reload before continuing on its route. Satellites like the one shown upper left in Fig. 1 are not used for synchronization purposes. Figure 2 shows the main challenge in our routing problem, namely the spatial and temporal synchronization of cargo bikes and vans. The abscissa in the diagram is the timeline and the ordinate shows the vehicles operating on a certain route. Van 1 starts at the van depot (large rectangle) and services two customers (dots). Then it arrives at the satellite (triangle) which represents a cloned satellites, where the cargo bike arrives after starting its route at the bike depot (small rectangle). In the meantime van 2 has also started its route from the van depot and services five customers before it arrives at the next satellite which represents also a cloned satellites and can therefore Table 1 Notation

Vehicle costs of a van/bike per unit of distance traveled
Distance from node i to node j for a van/bike Travel time of a van/bike from node i to node j Waiting time of a van/bike at node i Load of a van/bike after serving node i be located at the same physical location than the satellite mentioned above, but it can also be located at another physical satellite location. Here it meets with the cargo bike which has to be reloaded after servicing three customers.
To sum up the potential usage of a physical satellite s we can distinguish the following cases: 1. A physical satellite s is used once in the solution for the meeting of one bike with one van. Then exactly one of the respective cloned satellitess is used in the solution. 2. A physical satellite s is used once in the solution for the meeting of several bikes with one van. Then one respective cloned satellites for each bike is used in the solution. 3. A physical satellite s is used several times in the solution. Then for each usage in time either case 1. or case 2. can be realized. 4. A physical satellite s is not used at all in the solution, so all respective cloned satellitess stay unused as well.
Our problem is a generalization of the classical VRP and is thus NP-hard.

Solution framework
For our city distribution problem with spatial and temporal synchronization between vehicles from and delivery tours to customers on both echelons, we propose a solution procedure which basically consists of GRASP with path relinking (PR). This metaheuristic is chosen because we search for a fast and easy to use algorithm, which can help to decide, if such a sophisticated distribution scheme can be favorable in contrast to a typical one, where all goods are delivered by one type of vehicle to all customers inside and outside of the city center. Since GRASP is a greedy, memoryless procedure, we combine it with PR. This allows us to combine characteristics of good solutions as part of an intensification strategy. GRASP was first proposed by Feo and Resende (1989). In general, in each iteration a solution is generated by adding elements from a so-called restricted candidate list (RCL). Usually, the RCL is based on a greedy selection and an element from the list is selected randomly. Another way to combine greediness and random selection is to build the RCL of randomly selected elements and to choose the next element greedily based on costs. Once the solution is generated, it undergoes a local search phase.
The best solution found is stored during the search process and returned at the end (see Burke and Kendall (2005)). GRASP is often used in combination with PR to improve the performance of the algorithm (Resende and Ribeiro 2005). The initial solution is gradually transformed towards the guiding solution and the solutions on these paths are explored. In PR, the paths leading from one elite solution (ES) to a different one are explored as an intensification strategy. We maintain a solution pool that keeps the best and most diverse solutions. From this pool, two solutions, the initial and the guiding solution, are selected.
In the following subsections we describe our implementation of GRASP with PR.

Two-phase GRASP
For the construction of van routes on the first echelon we use information on demand and selection of potential satellites given by the bike routes on the second echelon. We thus have a two-phase GRASP, where we start with the construction of cargo bike routes.

Construction of cargo bike routes
The initial routes for cargo bikes are constructed by a nearest neighbor heuristic in which the next neighbor is chosen randomly from the RCL. RCL is a list of all possible (not routed) customers ordered by their insertion costs and is of size α. One customer of this list is either chosen randomly with uniform distribution (unbiased randomization), or randomly with geometric distribution with parameter β (biased randomization), i.e. the probability of being chosen decreases with increasing rank in the list (Juan et al. 2014).
Routes are constructed for each cargo bike successively. The starting node is chosen with (un)biased randomization from a list of satellites ordered by increasing proximity to all remaining bike customers. Then the nearest neighbor heuristic is performed as described above. The bike route is constructed until no more customers can be feasibly added with respect to the route duration constraint. The procedure is then performed for the next bike until all bike-customers are routed.
In order to preserve sufficient capacity and time in the route construction phase we use dummy satellites. This is an efficient way to avoid generating infeasible solutions. In cargo bike routes, we insert a dummy satellite whenever the cargo bike's capacity is exceeded. The cheapest satellite is used as a dummy satellite. At the end of the construction step, we remove all dummy satellites from the routes and perform a local search phase.
In the local search phase we use three operators: 2opt for intra-route changes, relocate for inter-route changes of two nodes and swap for intra-and inter-route changes of two nodes. This local search is performed for bike routes without any satellites, i.e., before the satellite insertion step.
The 2opt operator removes the connections between nodes (i, i + 1) and ( j, j + 1) and reconnects the two remaining parts of the route, i.e., the new connections are (i, j) and (i + 1, j + 1). As a consequence of this move the route segment between nodes j and i + 1 is reversed.
The relocate operator removes one node from one route and tries to insert it into another route. Note that this operator can remove the only customer of a route. This empty route is then deleted because it is no longer necessary.
The swap operator exchanges a node with another one either located in the same route or in a different route. All inter-route changes only consider routes of the same vehicle type and are based on best improvement. The operators described above are used as long as there is an improvement or the maximum number of local search iterations, θ , is reached.
After the local search phase we insert the satellites in the cargo bike routes with a Dynamic Programming (DP) approach based on the work by Beasley (1983). Based on a bike route consisting of the bike depot node and a number of customer nodes we calculate the insertion costs for the cheapest satellite nodes between each pair of nodes. Then we construct a directed graph that consists of the bike depot node as start and end node, and the insertion satellite nodes. Arcs in the graph represent cargo bike trips between two satellite visits. We only generate arcs that reflect the cargo bike's capacity. We can calculate the shortest path in the graph to split the giant tour in feasible cargo bike trips. Whenever we visit a node in this graph, this corresponds to a visit to a satellites in the original route. Figure 3 illustrates an example. The upper part shows the bike route starting at the bike depot (small rectangle) which is followed by a number of customers (dots; customer demand is given in brackets) and the bike depot as end node. The triangles in between represent the cheapest insertion satellitess between each pair of nodes (satellite number-for better orientation-and insertion cost are given below each satellite). The lower part of the figure represents the directed graph, where only feasible arcs with respect to the bike's capacity (given as 16) are depicted. The bold line represents the shortest path and indicates the insertion of satellites s1, s2 and s4 into the bike route to achieve a least costly feasible route.
After this insertion step it may happen that in a route segment between two satellites or the depot and a satellite a detour occurs (see Fig. 4). Therefore, we have added a segmentwise 2opt procedure-2optseg-to improve those segments.  Construction of van routes Once the cargo bike routes are constructed, we obtain the total demand for each inserted satellite by summing up the demand of all bikecustomers on the cargo bike trip originating from that satellite. Moreover, we can calculate the arrival time of the cargo bikes for the respective satellites and use this information for the construction of van routes.
The starting node for all van routes is the depot. The routes are constructed for each van by using the nearest neighbor approach described above. Once the capacity or route duration constraint is about to be exceeded, the van returns back to the depot. For van route construction we also use the concept of dummy satellites. In the initial construction step of a van route all satellitess used by cargo bikes are inserted in a van route without considering the temporal synchronization constraint. This is done to ensure feasible routes. Once all van routes are constructed, the satellites are removed and a local search is performed as described above.
After the construction and improvement steps of the van routes we have to insert all satellitess used by bikes in the van routes with respect to the temporal synchronization constraint. We use a simple best-fit approach, in which we start with the satellite with the earliest scheduled arrival time of a cargo bike. We search for the best (least costly) feasible insertion position in a van route.
At a satellites during the meeting between a bike and a van waiting times may occur and we can distinguish three cases: 1. The bike and the van arrive at the respective satellites exactly at the same point in time, then no waiting time is required (w b s = w ṽ s = 0). 2. The bike arrives earlier at satellites than the van, then the bike has to wait till the arrival of the van and the respective waiting time is assigned to satellites with w b s = t ṽ s − t b s ≤ w max . 3. The van arrives earlier at satellites than the bike, then the van has to wait till the arrival of the bike and the respective waiting time is assigned to satellites with w ṽ s = t b s − t ṽ s ≤ w max . Therefore, in our setting we have a kind of flexible time windows at the satellitess.
Before we can insert the next satellite the 2optseg is applied to the respective van route where the last satellite has been inserted and then we have to update all concerned bike and van route data with respect to arrival time, demand and overall route duration. Furthermore, if the bike has to wait at the satellite in question all arrival times have to be updated at the successive nodes in the same cargo bike route. If no feasible insertion position can be found for the satellite in question, then a new van route is created which is a mere depot-satellite-depot tour. This is repeated until all satellites are inserted in a van route and afterwards the feasibility of the solution is checked. Infeasible solutions are discarded.

PR-path relinking
The GRASP is enhanced with a PR phase. We use PR during GRASP as an intensification strategy and as post-optimization. For PR, we keep a pool of solutions, P, of size ξ . Whether or not a solution from the GRASP phase or after a PR step is accepted to the pool depends on the objective value and the diversity of the solution. A solution s can only be accepted if the objective function value is of a certain quality with respect to the best solution in the pool,s * , i.e., f (s) < f (s * )(1 + γ ), γ > 0. The diversity of a solution is measured by comparing the number of bikes used, the number of vans used, the number of synchronizations and the absolute difference of the number of satellite visits for each physical location of a satellite. In order to be accepted to the pool, a solution has to have a minimum difference, ν, to each solution in the pool. A solution that is better than the best solution in the pool is always accepted.
When a new solution is accepted, and the size of the pool becomes larger than the maximum pool size, ξ , another solution has to be removed. In that case, we remove the solution that is most similar to the new one and also that has a cost worse than the new one.
The procedure we use for path generation is based on Nguyen et al. (2012). We use PR on the second echelon, i.e., on the cargo bike routes only. First, the different trips are concatenated, starting with the cargo bike route with the earliest arrival time, to form a giant tour.
The path generation proceeds as follows: For giant tour A, which is a copy of the initiating solution, and giant tour B, which is the guiding solution and remains unchanged, we first make sure that they are of the same length and include the same nodes, i.e., if A has different satellite visits to B, these missing satellites are added and vice versa. Next, we look for a customer with a different position in A than in B. This customer is swapped in A with the node that is on its according position in B. Once there are no customers with different positions, the satellite occurrences are aligned by swaps and substitutions.
Using two example solutions with satellites s1−s3 and customers c1−c5 the PR procedure looks as follows.
Initial solution A: s1 c1 c2 c3 s2 c4 c5 Guiding solution B: s3 c4 c2 s2 c1 c3 s1 c5 As satellite s3 only occurs in B we have to add this satellite to A to balance the solutions. This leads to the balanced initial solution A: s1 c1 c2 c3 s2 c4 c5 s3 The PR starts with node c1 which is the first customer at position 1 in A, but its occurrence in B is at position 4. Therefore we swap c1 with s2. Proceeding this procedure leads to the following steps: s1 s2 c2 c3 c1 c4 c5 s3 s1 s2 c2 c4 c1 c3 c5 s3 s1 c4 c2 s2 c1 c3 c5 s3 s1 c4 c2 s2 c1 c3 s3 c5 Now all customers are swapped and we start with the satellites. As s1 is not at the correct position we need a swap and we get s3 c4 c2 s2 c1 c3 s1 c5 which is now identical with the guiding solution B.
After each swap we split the newly formed giant tour in single bike tours. Splitting is always done between a customer and a satellite before the maximum route duration is reached. Potentially occurring satellites at the end of a bike route are then discarded. If these bike routes are feasible, the above-mentioned GRASP steps to find a complete solution are added. That is, based on the information provided by the already constructed bike routes we generate then the van routes and insert the satellites as described in Sect. 4.1 with additional respect to the time synchronization constraint. For infeasible bike routes we insert new satellites to achieve feasibility and proceed as described. Feasible solutions are then checked for insertion in the ES pool.
PR is used each time a new GRASP solution is accepted to the non-empty pool. Then PR is performed between this solution and the solution from the pool which is most diverse. The solution with the better objective value acts as the initial solution. PR is also used as a post-optimization step when the pool size is reached. All pairs of solutions from the pool are combined and new solutions added to a new pool. If a new best solution is found, another iteration of post-optimization on this new pool is done and so on [see Resende and Ribeiro (2005)

Computational results
In this section, we present our computational results. First, we give a detailed instance description. Then we show the results for our parameter tests, and compare the performance of GRASP, GRASP with integrated PR, and GRASP with integrated PR and PR post-optimization phase. Finally, we present a comparison of the three different distribution policies described in Sect. 1.
The algorithm was implemented in C++ and tested under Linux Ubuntu 14.04 LTS running on a Virtual Machine (using two processors and 2 GB memory) on a host Intel(R) Core(TM) i5-3320 M CPU @ 2.60 GHz 4 GB RAM.

Instance description
For testing the performance of our algorithm we use three types of instances.
As first group of instances we use adapted Solomon instances based on those used by Grangier et al. (2014). The following additional changes are necessary:

Our van depot [resembling the CDC in Grangier et al. (2014)] is always located
at 50/100. 2. The depot in the Solomon instances is used as our bike depot. 3. We use the due time information of the Solomon-depot as our maximum route duration (only for the Solomon RC-instances this value has to be increased to enable feasible solutions). 4. As our two vehicle types drive with different speeds and have different costs we add these relations to the instance. We assume vans to drive three times as fast as cargo bikes. Furthermore vehicle costs for vans are set three times as high as for cargo bikes. The driver costs for van drivers are set 1.1 times higher than the costs for cargo bike drivers.
For more details about the primal adaptations we refer to Grangier et al. (2014). As the instances in each of the six groups (C1, C2, R1, R2, RC1 and RC2) only differ in the time windows, we use only the first instance out of each group (i.e. C101, C201, R101, R201, RC101, RC201) as test instances in our study. The second group of 12 instances is generated based on the idea of the Solomon instances applied to the city of Vienna. We use 2 clustered, 2 random and 2 combined settings with 100/125 customers and 10 satellites each. Each customer's demand d i is assumed as a homogeneous good and generated as a random integer between one and eight (i.e., half of the assumed maximum capacity of a cargo bike) as well as each customer's service time t l i , which is taken as a random integer between 6 and 16 min. We denote these instances by n100-c1, n100-c2, n100-r1, n100-r2, n100-rc1, n100-rc2, n125-c1, n125-c2, n125-r1, n125-r2, n125-rc1 and n125-rc2.
We also develop a real-world test instance based on the city of Vienna and therefore named 'vienna'. We choose the locations of 100 pharmacies at random and define the van depot as the location of one pharmacy wholesaler. The 18 satellites are placed randomly at potential locations (e.g. parking spaces) around the city center and the bike depot is located at a potential parking place for bikes in the city center (see Fig. 5).  Instead of euclidean distances we use travel times based on the real street network in Vienna. For the cargo bikes we use an average speed of 15 km/h (Energieregion 2009).
To calculate the travel times for the vans based on the street network we take Floating Car Data (FCD) from FLEET, an FCD system continuously working since 2003 and collecting data from about 3000 taxis in Vienna. Out of this data the system estimates the average speed per road element and time interval (Graser et al. 2012). We use the average speed-per-road element at 4:00 am on a typical weekday to calculate the travel times for the vans.
Service times and demands for each node are assumed as described above. The maximum route duration t max is assumed as 300 min for our test instances and as the duration of the time window of the depot for the Solomon instances. The maximum allowed waiting time w max is assumed as t max 10 which holds for all test instances. For the two homogeneous vehicle fleets we assume the values in Table 2 based on usual vehicle sizes and costs for vehicles and drivers. Furthermore, we use a penalty cost of 50 units to avoid van trips crossing the city center, which is defined as the largest possible in-circle C I with its center at the centroid of the polygon P V S formed by all satellites (see dashed circle in Fig. 5). These penalty cost is also used for the adapted Solomon instances described above.

Parameter settings
To define the parameters for our algorithms we perform tests for all test instances described in Sect. 5.1 with different parameter combinations (length of RCL α, maximum number of local search iterations θ , parameter of the geometric distribution β, size of the ES pool ξ , and diversity measure ν) and a fixed number of GRASP iterations η = 10. The cost threshold for accepting solutions to the ES pool γ is set to 1 in general to guarantee a sufficient number of diverse solutions.
Then we rank the results with respect to average total costs (column costs) in cost units and average computational time (column time) in seconds (see Table 3). Weighing both rankings with 80 % share of cost ranking and 20 % share of time ranking yields PC3/PC9/PC12 for GRASP, PC12/PC9/PC5 for GRASP+PR(I) and PC9/PC12/PC6 for GRASP+PR(I+P) as first/second/third place respectively. Based on these results we decide to use PC9 for all following tests in all settings (see Table 4).
For measuring the influence of runtime on the performance of each algorithm we perform tests with increasing number of iterations η. The maximum runtime of each test is limited with 300 s. The results as average values for all instances are shown in Fig. 6. This time limit of 300 s represents approximately η = 3000 for GRASP, η = 300 for GRASP+PR(I) and η = 190 for GRASP+PR(I+P). Based on the convergence behavior of our algorithms we decide to use a maximum of η = 500, 200, 100 for GRASP, GRASP+PR(I) and GRASP+PR(I+P) respectively and limit the runtime additionally at 300 s.

Results for testing GRASP and PR components
We have solved all instances using -GRASP alone, -GRASP+PR(I), where the PR step is used as integrated step, and -GRASP+PR(I+P), where the PR step is also used as a post-optimization phase.
All instances have been solved five times and the average value for each instance is presented in all tables. Table 5 depicts the aggregated total costs for the all instances. Values in brackets show the improvement of each method compared to the previous column.
Detailed results for all instances can be found in Tables 6, 7 and 8. Here the number of used bikes/vans (#b/v), the distance traveled by all bikes/vans (dist b/v), the occurring waiting time for bikes/vans (wait b/v), costs for bikes/vans (cost b/v), the sum of all costs ( cost), as well as the number of used synchronization points (#sync) is depicted.  Our computational results illustrate that for most of the instances with randomly assigned customers the PR post-optimization step does not lead to an improvement at all, while in general, instances with clustered customers can be further improved by this step. Besides, the more synchronizations occur in each bike tour, the better the PR step performs which can be seen for example at the results of instance Vienna, where 2 bikes are used on average and a total of 5 synchronizations is necessary. Here the PR(I) step can improve the result by about 3 % and the post-optimization step can add about another 7 % improvement to the result (see Tables 6, 7 and 8 as well as Table 5). If each bike tour requires only one synchronization in a solution, the improvement gained by the PR step is on average rather small (e.g. instances R201 and RC201).

Comparison of different distribution policies
We investigate three different distribution policies. The first one resembles the current distribution policy classically used by companies which is the solution of a capacitated VRP. In the second one, the use of intermediate satellite facilities is possible without additional costs. These satellite facilities are assumed to have storage space, so that a temporal synchronization between the first and second echelons is not necessary. Finally, the last distribution policy is the two-echelon distribution scheme with temporal synchronization of vans and cargo bikes investigated in detail in this paper. As our PR builds on bike routes, we use only GRASP to solve the first distribution policy, where only vans are used. The other two distribution policies are solved by GRASP+PR(I+P). Table 9 presents the aggregated results for the three distribution policies, where numbers in brackets refer again to the increase (+) or decrease (−) in costs as a percentage value compared to the previous column. Please note that we do not consider storage cost for the second distribution policy. Detailed results can be found in Tables 10, 11 and 12 depicting the same columns as in Tables 6, 7 and 8. Table 9 shows that in instance R201 the combined usage of cargo bikes and vans (column labeled 'without sync') is less costly than the usage of vans alone, although for the solution with vans alone we do not even use a penalty cost for van trips crossing the city center. This can be explained by the fact that in this instance one van can be substituted by two bikes (see Tables 10 and 11) and the underlying cost structure (see Table 2).
For all other instances the combined usage of cargo bikes and vans causes higher costs, but since some vans can be substituted by cargo bikes, a significant reduction of CO 2 -emissions can be provided and the presence of vans in the city center can be reduced. Policies with temporal synchronization (column labeled 'with sync' in Table 9) always lead to a cost increase compared to policies without synchronization. The level of cost increase for the synchronization step depends mainly on the number of required temporal synchronizations compared to the number of required bikes (see Tables 10,  11 and 12). These synchronizations can lead to additional waiting times of vans and/or cargo bikes if the required synchronizations do not only occur at the beginning of a bike route but also in between (i.e. a bike meets with a van not only after leaving the bike depot but also afterwards as in instance n125-c1, where about 3 bikes require about 9 synchronizations). Besides, it can also happen that an additional vehicle is required to achieve a feasible solution as it can be seen for example at instance n100-rc1, where instead of 8 vans now 9 vans are necessary to achieve a feasible solution in case of temporal synchronization. However, given that we do not consider the costs of storage facilities in the 'without sync' policy of Table 9, the additional cost of synchronization are quite low on average with 3.15 %, although those values differ depending on the instance in the range of 0.23 to about 10 %. Comparing the additional costs for synchronization to the storage costs can help planners decide which distribution policy is more preferable. Figure 7 shows a sample solution for our Vienna test instance where six vans are used to deliver all occurring demand assuming that there are no restricted zones Table 10 Detailed results for distribution policy with vans only  and no additional costs in the city center, whereas Fig. 8 displays a sample solution where six vans (solid lines) and two cargo bikes (dashed lines) are used and temporal synchronization takes place at satellites (triangles). Here we penalize van trips crossing the city center as mentioned before.

Conclusion
In this paper we present an innovative city distribution scheme for a two-echelon routing problem with temporal and spatial synchronization between vans and cargo bikes. Van trips in the city center should be avoided to mitigate the impact of traffic congestion and noise nuisance. Therefore, cargo bikes perform the last mile distribution in the city center.
We present a GRASP with PR as solution method. Our algorithm provides a first and fast support for decision makers to evaluate when it is-even economicallypreferable to use green modes of transport in a city distribution network. Especially the low computational costs of solving this complex problem have to be mentioned, as shown by running our algorithm on a Virtual Machine. We test different possibilities for the inclusion of PR in GRASP. Computational results show that GRASP with integrated PR outperforms GRASP alone in all instances. Using PR also as a post-optimization phase leads to slightly better results in nearly all cases, though the computational time increases significantly. A comparison of the three distribution policies illustrates that for some instances costs can be saved by the combined usage of cargo bikes and vans instead of vans alone. In all cases emissions can be reduced through the substitution of vans by cargo bikes. In all cases the policy which includes synchronization is more costly than the one without synchronization, when storage costs are not taken into account. Planners can evaluate the trade-off between the additional operational costs of synchronization and the storage costs of satellite facilities.
Nevertheless, certain technical requirements for cargo bikes such as ensuring the required temperature range in the cargo may have to be tackled to use our distribution scheme in a real-world application.
Future research will include the evaluation of C O 2 -emissions as well as uncertainties in travel times and dynamic requests.