A variable neighborhood search algorithm with constraint relaxation for the two-echelon vehicle routing problem with simultaneous delivery and pickup demands

This paper considers a special vehicle routing problem, the two-echelon vehicle routing problem with simultaneous delivery and pickup demands (2E-VRPSDP). The 2EVRPSDP differs from classic transportation and vehicle routing problems in two ways. First, freight delivery from the depot to the customers is managed by shipping the freight through intermediate satellites. Second, each customer in the 2E-VRPSDP may have simultaneous delivery and pickup demands. The 2E-VRPSDP is an extension of the two-echelon vehicle routing problem (2E-VRP) and the vehicle routing problem with simultaneous delivery and pickup (VRPSDP). A variable neighborhood search algorithm is designed to solve the 2E-VRPSDP in which both feasible and infeasible solutions can be explored. Numerical results show that the proposed algorithm is effective and that the algorithm can provide reasonable solutions within an acceptable computational time.


Introduction
In the past decade, the two-echelon vehicle routing problem (2E-VRP) has become a new interest in the vehicle routing problems (VRPs) research field. The logistics network of 2E-VRP is composed of two echelons, i.e., the freight is delivered from the depot to the satellites by large capacity vehicles in the first echelon and then transported from the satellites to the customers by the second echelon vehicles of relatively small capacity. The objective of the 2E-VRP is to minimize the total transportation costs of both echelons. The 2E-VRP problem arises in many practical transportation and distribution contexts, such as city logistics applications (Feliu et al. 2007;Jepsen et al. 2013).
The 2E-VRP has been studied by many researchers, and some very good studies can be found for this problem Jepsen et al. 2013;Breunig et al. 2016;Liu et al. 2017). In the existing research, generally, it is assumed that one customer only has either delivery or pickup demand. However, in some practical applications, the vehicle needs to distribute and collect freight simultaneously. We consider a practical two-echelon logistical problem arising in the home healthcare (HHC) industry in Shanghai, China. Many HHC companies have been built in Shanghai for patients who require long and regular health care to provide quality health services at their homes. The typical services in the HHC involve various logistic activities, including delivering medicines and medical instruments to patients, picking up biological samples from patients' homes, and collecting medical waste from patients' homes for disposal. Considering the HHC company, the core component in-home healthcare problems is to find a feasible working schedule for their drivers and vehicles to reduce the operating cost. For a large city, because customers spread in a very large area, one HHC company establishes some subcompanies in different districts, and daily logistic services such as delivering medicines and collecting medical waste are implemented as a two-echelon logistics system. First, in the morning, the goods to be delivered to the customers (e.g., the medicines for the customers) are transported from the company (or center warehouse) to the subcompanies. Next, small subcompany vehicles deliver goods to customers and pick up goods from customers (e.g., collecting medical waste) to take back to the subcompanies. After all the pickup goods are collected in subcompanies, finally, the vehicles at the HHC company pick up such goods from subcompanies to return to the HHC company and complete one day's operations. Clearly, this is a special two-echelon transportation network in which the first echelon consists of the HHC company and subcompanies, and the second echelon contains subcompanies and customers. We can further observe that, unlike the classical 2E-VRP, the operations by the HHC company consist of three phases. The HHC company must design a set of routes for vehicles in this two-echelon network. According to the HHC company investigations, it is not easy for the company's manager and planner. This is because this routing optimization problem is rather complex, which can be seen as a combination of two challenging problems, i.e., the two-echelon vehicle routing problem (2E-VRP) and the vehicle routing problem with simultaneous delivery and pickup (VRPSDP), because the problem has a two-echelon network structure and the customers have both delivery and pickup demands simultaneously. Since both the VRPSDP and the 2E-VRP are NP-hard, the problem stated above is also NP-hard and even more complex and harder than these two problems.
In this paper, we focus on this special logistical problem, which is named the two-echelon vehicle routing problem with simultaneous delivery and pickup (2E-VRPSDP). Although the 2E-VRP and VRPSDP have been studied for more than ten years, as stated above, the 2E-VRPSDP solution structure differs from the classical problems (the detailed differences are shown and analyzed in Sect. 2.) To our knowledge, there is no literature about the 2E-VRPSDP proposed in this paper. Because of the difficulties in solving instances of practical interest, we propose a variable neighborhood search (VNS) heuristic in which both feasible and infeasible solutions are explored to enhance the algorithmic search ability. Numerical results show that the VNS algorithm can efficiently solve the 2E-VRPSDP.
The rest of this paper is organized as follows. Section 2 defines the problem and discusses the differences with existing problems that have been studied in the literature. Section 3 builds the mathematical model for the 2E-VRPSDP. Section 4 reviews the relevant literature. Section 5 describes the VNS algorithm. Computational results are reported in Sect. 6. Section 7 concludes the paper.

Problem definition
The 2E-VRPSDP can be defined formally as follows. Let G = (V, E) be a graph with the node set V = V 0 [ V S [ V C and edge set E = E 1 [ E 2 . In set V, V 0 = {0} is the depot, set V s = {1, …, |V S |} represents the set of satellites where |V S | is the cardinality of set V S , and set V C = {|V S |? 1, …, |V S | ?|V C |} is the set of customers. Each customer i [ V C has a delivery demand and a pickup demand simultaneously, denoted as d i and p i , respectively. A fleet of K 1 homogeneous first echelon vehicles with a capacity of Q 1 is located at depot V 0 , which can visit the depot and satellites. Additionally, a total of K 2 homogeneous second echelon vehicles with a capacity of Q 2 are located at all satellites, which can only visit the satellite and customers, and Q 1 [ Q 2 . Set E is divided into two subsets, E 1 and E 2 . Set corresponds to the edges between the depot and the satellites and between different satellites, which can only be traveled by the first echelon vehicles. Set E 2 ¼ f i; j ð Þ : i; j 2 V S [ V C ; i; j 6 2 V S Â V S g represents the edges connecting the satellites and the customers and those connecting different customers. Each edge (i, j) [ E has a nonnegative cost (distance), denoted as c ij , and for each pair Note that a solution of the 2E-VRPSDP consists of three phases. In the first phase, a number of first echelon vehicles start from depot V 0 and deliver goods to a sequence of satellites and then return to the depot. For notational convenience, each of these routes is called the first phase delivery route. Then, in the second phase, second echelon vehicles start from the satellites, each meeting the delivery and pickup demands of one or several customers, and return to the starting satellites. Now, all the pickup demands from the customers are stored at the satellites. Similarly, each route at this echelon is referred to as a second phase delivery and pickup route. Finally, again, the first echelon vehicles depart from the depot and collect all pickup goods from the satellites and return back to the central depot. Such routes are called third phase pickup routes. Note that at each phase, direct vehicle routes between the depot and customers are forbidden. At each phase, some constraints must be respected as follows.
(1) Each first phase delivery route and third phase pickup route start and end at the depot; (2) Each route is assigned to exactly one vehicle; (3) At any time, the total load of the vehicle cannot exceed the vehicle capacity; (4) A customer can be served by only one second phase delivery and pickup route; a satellite can be served by at most one first phase delivery route and at most one-third phase pickup route. This constraint means that the demands of both echelons cannot be split.
The objective of 2E-VRPSDP is to determine the vehicle routes of both echelons to the sum of the routing costs, i.e., the total travel distances of all vehicle routes in the above three phases are minimized.
Compared with the 2E-VRPSDP, the 2E-VRP has only two-phase operations: first, the goods to be delivered to the customers are transported from the depot to the satellites by large vehicles; next, small vehicles starting from satellites deliver goods to customers and complete the operations (Feliu et al. 2007;Jepsen et al. 2013;Hemmelmayr et al. 2012). To clearly illustrate the 2E-VRPSDP and compare it with the 2E-VRP, a solution to the 2E-VRP with one depot, three satellites (S 1 -S 3 ), and seven customers (C 1 -C 7 ) is shown in Fig. 1, which consists of three parts A-C. Each customer has pickup and delivery demands. First, as shown in part A of Fig. 1, the company distributes the goods to the satellites by two first echelon delivery routes (blue lines). Next, in part B, three vehicles start from each satellite and deliver goods to each customer and collect pickup goods to deliver back to the satellites (black dotted line). Finally, part C shows two routes that collect all the goods to deliver to the depot (red lines). Belgin et al. (2018) introduced an interesting two-echelon vehicle routing problem with simultaneous pickup and delivery requests. In their problem, the customers have pickup and delivery demands, and the pickup and delivery activities are performed simultaneously by the same vehicles through satellites to customers in the second echelon. Belgin et al. (2018) assumed that in the first echelon, goods are delivered from the central depot to each satellite and collected from satellites to deliver back to the depot by the same vehicles. That is, a solution to the problem of Belgin et al. (2018) consists of only two steps: The delivery routes and pickup routes between the depot and satellites (as shown in parts A and C of Fig. 1) are combined. Furthermore, we note that many other researchers also state that pickup and delivery in two-echelon city logistics can improve the efficiency if it is performed simultaneously on the first echelon vehicles. For example, Gianessi et al. (2016) solved the related tactical planning problems; for tactical planning problems, please refer to Fontaine et al. (2017Fontaine et al. ( , 2021. We agree that the operations in these works are reasonable in practical applications. We investigate many logistics and home healthcare companies and find that they do not adopt a two-step model because, in some scenarios, pickup and delivery operations for a satellite are not easy to perform at the same time. For example, one first echelon vehicle a reaches satellite i at 8:00 and unloads the goods for this satellite. Next, a second echelon vehicle b delivers such goods to customers and returns to this satellite at 10:00, with all pickup demands from its customers. Following the assumption of Belgin et al. (2018), vehicle a has to wait at satellite i until vehicle b returns to this satellite, i.e., two hours in this example. In practice, some companies deliver goods from a central depot to satellites in the morning and collect goods from satellites to deliver back to the depot in the afternoon, consistent with the problem definition in this paper.

Mathematical formulation
The mathematical formulation of the 2E-VRPSDP (mixedinteger linear program, MILP) is presented below.
Decision variables: y i,j binary variable equal to 1 if a vehicle travels directly from node i to node j at a first phase delivery route; x i,j binary variable equal to 1 if a vehicle travels directly from node i to node j at a second phase delivery and pickup route; l i,j binary variable equal to 1 if a vehicle travels directly from node i to node j at a third phase pickup route; z i,j binary variable equal to 1 if customer i is assigned to satellite j; a i quantity of goods to be delivered and loaded on the second phase vehicle until customer i is visited; b i quantity of pickup goods loaded on the second phase vehicle immediately after customer i is visited; n i quantity of demands to be delivered at satellite i; o i quantity of demands to be picked up at satellite i; A i quantity of goods to be delivered to satellites loaded on the first phase vehicle until satellite i is visited; B i quantity of pickup goods loaded on the third phase vehicle immediately after satellite i is visited; Objective function: x i;j z i;j 8i 2 V C ; j 2 V S ð10Þ x e;m þ z e;i þ X j2V S ;j6 ¼i z m;j 2 8e; m 2 V C ; e 6 ¼ m; The objective function (1) minimizes the total transportation cost. Constraints (2) ensure that each customer is visited exactly once. Constraints (3) guarantee flow equations for customers, i.e., every vehicle that arrives at a customer must leave that customer. Vehicle load constraints in the second phase delivery and pickup route are explained in Constraints (4)-(8). Such constraints also eliminate the subtours in the second echelon routes. Constraint (6) ensures that in any customer's area after the vehicle delivers and picks up goods, the loaded quantity does not exceed the vehicle capacity. Constraint (7) denotes that the demand of customer i must not exceed the load quantity on the vehicle before arriving, and this load quantity value must not exceed the capacity of the vehicle.
In Constraint (4) if x i,j = 1, i.e., customer j is next to customer i, and a id i = a j , which indicates the delivery quantity on the second phase vehicle is decreased by the demand quantity of customer i after visiting, so that a j þ d i þ Q 2 Q 2 þ a i . Otherwise, when x i,j = 0, the inequality also holds since d i B a i and a j B Q 2 according to the definition of these notations. Constraints (5) and (8) correspond to the pickup goods case and are of a similar form to Constraints (4) and (7). Constraints (9) force each customer to be assigned to one satellite. Constraints (10)-(11) ensure that there is no arc connecting customer i and satellite j if customer i is not allocated to satellite j. Constraint (12) avoids two connected customers belonging to two different satellites. Constraints (13) and (14)

Literature review
Although many variants of the VRP and the 2E-VRP have been studied in the existing literature, little research has been performed on the 2E-VRPSDP. As stated above, two main bodies, the 2E-VRP and the VRPSDP, are relevant to the 2E-VRPSDP. We survey the literature in these two bodies of research.

The 2E-VRP literature
The 2E-VRP has drawn much attention in the past decade. Both exact and heuristic algorithms have been designed to solve this problem. In terms of exact algorithms, first, index formulation similar to that of Perboli et al. (2011); they derived a relaxation from it to avoid giving incorrect upper bounds when more than two satellites were included in the solution. The branch-and-cut algorithm of Jepsen et al. (2013) performed better than that of Perboli et al. (2011). A branch-and-cut-and-price algorithm was presented by Santos et al. (2014), which overcame symmetry issues and was further strengthened by valid inequalities. Baldacci et al. (2013) proposed a new mathematical formulation. Their exact method decomposed the problem into a limited set of multiple-depot vehicle routing problems with side constraints. Computational results on benchmarks proved that their method outperformed previously published exact methods in terms of size, number of problems solved to optimality, and computing time. Considering some variants of the classical 2E-VRP, Liu et al. (2018) introduced a 2E-VRP with grouping constraints, in which customers were divided into several disjoint groups, and the grouping constraints ensured that customers from the same group were served by vehicles from the same satellite. They formulated the problem as a mixed-integer program and proposed valid inequalities to strengthen the model. A branch-and-cut algorithm was implemented to solve the problem, which could solve more instances to optimality than CPLEX. Darvish et al. (2019) studied a flexible 2E-VRP in which a supplier delivered a commodity to its customers through a two-echelon supply network. Two sources of flexibility were analyzed: flexibility in network design and flexibility in due dates. The former was related to the possibility of renting any of the satellites in any period of the planning horizon, whereas the latter was related to the possibility of serving a customer between the period an order was set and a due date. A mathematical programming formulation to this problem was presented, and an exact method was proposed that was based on the interplay between two branch-and-bound algorithms. Dellaert et al. (2019) studied the 2E-VRP with time windows. Different from the classical 2E-VRP, the second echelon consisted of transferring freight from satellites to the final customers within their time windows. They proposed two path-based mathematical formulations for the problem and developed branch-and-price-based algorithms to solve the problem. The algorithms solved instances of up to five satellites and 100 customers to optimality. Breunig et al. (2019) formulated an extension of the 2E-VRP, called the electric 2E-VRP, which involved electric vehicles for second echelon deliveries, battery capacity constraints, and possible visits to charging stations, and used it as a prototypical problem for the study of multiechelon batterypowered supply chains. They designed an efficient exact algorithm based on the enumeration of candidate solutions for the first echelon and on bounding functions and route enumeration for the second echelon, along with a problem-tailored large neighborhood search metaheuristic. Marques et al. (2020) proposed a branch-cut-and-price algorithm for the 2E-VRP. The authors introduced a new route-based formulation for the problems that do not use variables to determine product flows in satellites. They also introduced a new branching strategy that significantly decreased the size of the branch-and-bound tree and introduced a new family of satellite supply inequalities. Because exact algorithms usually cannot effectively solve the large-size 2E-VRP in a reasonable computation time, many researchers have resorted to heuristic methods. Crainic et al. (2011) developed a multistart heuristic based on separating the problem by solving customer assignments heuristically and then dealing with the remaining VRPs. In their method, a perturbation mechanism was adopted to iteratively build new solutions, and a feasibility search procedure was used to bring the solution back into the feasible region. In addition to the exact methods, Perboli et al. (2011) also presented two math heuristics based on the information obtained by solving the linear relaxation of the mathematical formulation model. Hemmelmayr et al. (2012) designed an adaptive large neighborhood search heuristic for the problem. They introduced some largescale instances with up to 200 customers. Zeng et al. (2014) solved the problem using a greedy randomized adaptive search procedure embedded with a route-first cluster-second procedure and a variable neighborhood descent. Their approach was tested on instances not larger than 50 customers. Breunig et al. (2016) designed an effective large neighborhood-based heuristic for solving the 2E-VRP. Their algorithm improved 18 best-known solutions. Grangier et al. (2016) addressed a variant of the 2E-VRP that integrated constraints arising in city logistics such as time window constraints, synchronization constraints, and multiple trips at the second echelon. They proposed an adaptive large neighborhood search to solve this problem. Amarouche et al. (2018) proposed a new hybrid heuristic method for solving the 2E-VRP that relied on two components. The first component effectively explored the search space to discover a set of interesting routes, and the second recombined the discovered routes into high-quality solutions. Rohmer et al. (2019) presented a two-echelon inventory routing problem for perishable products. Products were delivered from a supplier to an intermediary depot, where storage may occur and from which they were delivered by smaller vehicles to the customer locations.
The objective was to minimize the total transportation and holding costs. An adaptive large neighborhood search metaheuristic was designed to address the problem. Jie et al. (2019) considered the 2E-VRP with battery swapping stations, which aimed to determine the delivery strategy under battery driving range limitations for city logistics. The electric vehicles operating in the different echelons had different load capacities, battery driving ranges, power consumption rates, and battery swapping costs. The authors proposed a hybrid algorithm that combined column generation and an adaptive large neighborhood search to solve the problem. Mühlbauer and Fontaine (2021) studied a variant of the 2E-VRP that used cross-docking from vans to cargo bicycles at so-called satellites. To solve large instances, the authors introduced a new efficient parallelized large neighborhood search algorithm. The algorithm was tested using symmetric 2E-VRP benchmark instances from the literature.

The VRPSDP literature
In addition to the 2E-VRP, the VRPSDP has also been studied intensively (Koç et al. 2020). The VRPSDP was first proposed by Min (1989 (2021) addressed a generalization of the one-commodity pickup and delivery traveling salesman problem where each customer supplied or demanded a given quantity of a certain product. The authors presented three mathematical models for the problem and designed an exact branch-and-cut algorithm to solve it. Compared with the exact algorithms, heuristic approaches dominated the solution methodologies of the VRPSDP. Nagy and Salhi (2005) proposed a heuristic that allowed infeasibilities to occur and guided the search toward strong feasibility through search routines. Crispim and Brandao (2005) presented a hybrid algorithm for the VRPSDP that was comprised of the two metaheuristics of tabu search and variable neighborhood descent. Alfredo Tang Montané and Galvão (2006) proposed a tabu search approach that utilized relocation, interchange and crossover movements to obtain interroute adjacent solutions and used a 2-opt procedure to obtain alternative intraroute solutions. Bianchessi and Righini (2007) presented and compared constructive algorithms, local search algorithms and tabu search algorithms on the VRPSDP. Gajpal and Abad (2009) designed an ant colony system (ACS) for solving the VRPSDP, which used a construction rule as well as two multiroute local search schemes. Zachariadis et al. (2010) introduced an adaptive memory algorithmic framework to solve the VRPSDP, which combined promising solution features to generate high-quality solutions. Avci and Topaloglu (2015) proposed an adaptive local search solution approach for the VRPSDP, which hybridized a simulated annealing inspired algorithm with variable neighborhood descent. A perturbation-based variable neighborhood search heuristic for solving the VRPSDP was designed by Polat et al. (2015). Zachariadis et al. (2016) introduced a special VRPSDP with two-dimensional loading constraints, which covered cases where customers posed delivery and pick up requests for transporting nonstackable rectangular items. The authors proposed an optimization framework that employed memorization techniques to accelerate the solution methodology. Computational results were reported on the VRPSPD, the VRP with two-dimensional constraints, and newly constructed benchmark problems. Qiu et al. (2018) studied a variant of VRPSDP, in which customers' demands were discrete in terms of batches (or orders) and each customer could be visited by a variety of vehicles or several times by one vehicle. A tabu search algorithm with specifically designed batch combinations and item creation operations was proposed to solve the problem. Majidi et al. (2018) dealt with the pollution routing version of VRPSDP, where the goal was to minimize fuel consumption and emissions by scheduling and routing customers. A nonlinear mix integer programming model was presented for this problem, and an adaptive large neighborhood search heuristic was proposed for the solution method including new removal and insertion operators. Zhang et al. (2019) addressed a multicommodity many-to-many vehicle routing problem with simultaneous pickup and delivery for a fast-fashion retailer in Singapore. To solve large-scale instances, an adaptive memory programming-based algorithm combined with techniques such as the regret insertion method for initializing the solution pool, the segment-based evaluation scheme, and the advanced pool management method was proposed in this work.
Although as stated above, many good exact algorithms and heuristics were proposed for the 2E-VRP and VRPSDP, to our knowledge, only Belgin et al. (2018) studied a special two-echelon vehicle routing problem with simultaneous pickup and delivery. However, the basic operation scheme of the problem in Belgin et al. (2018) differs from the problem studied in this paper. Belgin et al. (2018) assumed that goods of the first echelon are delivered from the central depot to each satellite and collected from satellites back to the depot simultaneously by the same vehicles. As stated in Sect. 2, our 2E-VRPSDP separates the routes in the first echelon into one delivery route and one pickup route. To our knowledge, this is the first work on 2E-VRPSDP.

Variable neighborhood search algorithms
Since the VRPSDP and the 2E-VRP are two difficult NPhard problems, most research on these two problems resorts to heuristics methods, such as tabu search (TS), variable neighborhood search (VNS), adaptive large neighborhood search heuristic (ALNS), and genetic algorithm (GA). Because the 2E-VRP can be seen as a special case of the 2E-VRPSDP studied in this paper, the 2E-VRPSDP is also an NP-hard problem and even more complex than VRPSDP and 2E-VRP. Thus, exact algorithm or commercial solver behavior is highly unpredictable. For example, we attempt to solve the above 2E-VRPSDP mathematical formulation using the Cplex and Gurobi. We find that the solvers cannot find the high-quality solution in a reasonable computational time (e.g., 2-3 h) even for moderate-size test instances (refer to the numerical results in Sect. 6.4). Thus, in this paper, a variable neighborhood search algorithm (VNS) is designed for the 2E-VRPSDP. The principle of VNS is a systematic change in neighborhoods within a local search procedure. As a local search-based algorithm, the main advantage of VNS is not the computation time. For example, it is slower than some simple greedy or hill-climbing algorithms. However, it has a great ability to find high-quality solutions, especially for some very complex NP-hard optimization problems. This has been widely proven by solving the VRP as well as its variants (Kytöjoki et al. 2007;Fleszar et al. 2009;Hemmelmayr et al. 2009). For a more thorough description of VNS, refer to Mladenović and Hansen (1997) and Hansen and Mladenović (2001). In this paper, we build a new VNS approach, the outline of which is shown in Algorithm 1. First, the algorithm identifies a set of neighborhood structures N k (k = 1, 2, …, k max ) and an initial solution S. Parameter k equals 1. Then, a shaking step is performed by randomly selecting a solution S 0 of solution S based on the kth neighborhood N k . Next, the algorithm applies a local search improvement procedure to the current solution S 0 . This procedure is repeated once a new incumbent solution S 00 is found. If this local optimum S 00 is better than S, the algorithm sets S / S 00 and searches with k = 1. Otherwise, VNS switches to the next neighborhood with k = k ? 1 to perform another shaking step, followed by local search improvement. Note that our approach extends the search to infeasible solutions, i.e., both feasible and infeasible solutions are allowed by the VNS. Each solution to the problem corresponds to a set of routes. Infeasibility occurs for a VNS solution if the total demand of a vehicle exceeds the vehicle's specified capacity. The algorithm uses a weighted linear penalty function for violations of this constraint. The details of the implementation of an infeasible solution are presented in the following sections

General functioning of the algorithm
Since the 2E-VRPSDP is a rather complicated problem, the VNS approach extends the search space to the infeasible region. It is frequently conjectured that infeasible solutions enable a better transition during the search between structurally different feasible solutions (Cordeau et al. 2001). In particular, purposeful management of penalties may enable us to focus the search toward borders of feasibility, a place where high-quality solutions are more likely to be located (Glover and Hao, 2011). Thus, in our VNS, the vehicle capacity constraint is temporarily violated during the search. A weighted and linear penalty function with a penalty weight a is introduced into the cost function for the violation of vehicle capacity constraints. For a solution S, a general cost function c(S) is defined as c(S) = d(S) ? a 9 e(S) for solution evaluation, where d(S) denotes the total travel distances of vehicle routes (original objective function), e(S) denotes the total violations of vehicle load of all the routes in three phases, and a is the penalization parameter. Term e(S) is defined as follows: where x ? = max(0, x), load i1 and load i3 are the vehicle loads after a first phase delivery route and a third phase pickup route serves satellite i, respectively; load i2 is the load of a vehicle after a second phase delivery and pickup route serves customer i. Clearly, if solution S is feasible, i.e., the vehicle capacity is strictly satisfied at each node, e(S) = 0 and c(S) = d(S); otherwise, e(S) [ 0 and c(S) [ d(S). In our VNS penalty parameter, a is adjusted dynamically to facilitate the exploration of the search space. If the solution S 00 obtained after shaking and local search is feasible (referred to Algorithm 1), the value of a is divided by a factor of 1 ? u (u [ 0); otherwise, a is multiplied by this factor. This penalization parameter a is initialized with a 0 and is limited between an upper bound a max and a lower bound a min , which limit the maximum and minimum values of this parameter during the search process.

Construction of an initial solution
Based on the savings algorithm (Clark and Wright 1964), we propose a method to obtain an initial solution to our problem. The initial solution method has three stages.
• Stage 1: customer assignment First, a random sequence of all the customers is generated. Following this sequence, each customer is assigned to its nearest satellite. During this assignment procedure, the algorithm ensures that for each satellite, the total delivery or pickup demands from all associated customers cannot be greater than the first echelon vehicle capacity Q 1 . If an assignment of a customer to a satellite violates this constraint, this customer is assigned to its second nearest satellite, and so on until all the customers are assigned to satellites. Because the vehicle capacity is limited, the above procedure may fail to assign each customer to a satellite. In this case, the initial customer sequence is regenerated until all customers have been assigned. Note that in the first step, the random sequence introduces randomness to the algorithm, similar to the works of Cordeau et al. (1997Cordeau et al. ( , 2001). • Stage 2: second echelon routing The algorithm applies the savings algorithm (Clark and Wright 1964) to each satellite to service the customers who have been assigned to this satellite. For satellite s, a separate back and forth route for each customer is initially created. For each customer pair i and j of satellite s, the savings is defined as s ij = c si ? c jsc ij . Starting from the largest and nonnegative saving s ij , two routes (s, …, i, s) and (s, j, …, s) are merged into a single second phase delivery and pickup route (s, …, i, j, …, s). The above route combinations are executed until the number of second phase delivery and pickup routes is no more than K 2 (the total number of second echelon vehicles). Note that the above procedures may construct an infeasible solution that violates the vehicle capacity constraint. For example, part A of Fig. 2 shows four routes that start and end at satellites S 1 and S 2 . Based on such routes, the saving values of the merging route (S 1 , C 1 , C 2 , S 1 ) and (S 1 , C 3 , S 1 ) and combining (S 2 , C 4 , S 2 ) and (S 2 , C 5 , S 2 ) are computed. If the former savings value is larger, then the two routes of satellite S 1 are merged, and one larger route (S 1 , C 1 , C 2 , C 3 , S 1 ) is obtained, as shown in part B of Fig. 2. • Stage 3: first echelon routing After obtaining the second phase delivery and pickup routes, similar to stage 2, the algorithms apply the savings algorithm to obtain the first phase delivery routes and the third phase pickup routes.

Shaking
Shaking is the foundation of the VNS, which diversifies the solution and avoids it from being trapped in a local optimum. Each neighborhood should determine a proper balance between perturbing the incumbent solution and retaining the good parts of the incumbent solution. The algorithm employs ten interroute neighborhood structures (namely, N 1 -N 10 ) in VNS shaking. The set of neighborhoods is summarized in Table 1. Neighborhoods N 1 to N 4 relocate or exchange satellite(s) among two first phase delivery routes. Neighborhoods N 5 to N 8 similarly deal with the third phase pickup routes. Note that if a satellite is moved, the customers and second phase routes associated with this satellite are moved with this satellite. Figure 3 shows an example of a neighborhood N 1 operator: Satellite S 1 is removed from a first phase delivery route (Depot-S 2 -S 1 -Depot) and inserted into another route (Depot-S 3 -Depot).
The last two neighborhoods, N 9 and N 10 , involve the change in the second phase delivery and pickup route and the first echelon routes (pickup route and delivery route), simultaneously. They are more complex than the above neighborhoods. To clarify N 9 and N 10 , two satellite statuses close and open are used. In this paper, when one or more customers are assigned to a satellite, the status of this satellite is open. For a satellite, when all its customers are assigned to other satellite(s), this satellite becomes closed. Neighborhoods N 9 first randomly close a satellite. The customers assigned to this satellite are removed from this satellite one by one. The algorithm inserts each removed customer into a random position of another second phase delivery and pickup route or constructs a new second phase delivery and pickup route that only contains this customer. In terms of constructing a new second phase route, it may be built on an open satellite or a currently closed satellite. Once this new route is built on a closed satellite, this satellite changes to open. Thus, this satellite is simultaneously inserted into a first phase delivery route and third phase pickup route. For example, in Fig. 4, after shaking, satellite S 2 is closed. Customers C 4 and C 5 are assigned to the current open satellite S 1 and the closed satellite S 4 , respectively. Thus, satellite S 4 now becomes open and is inserted into a first phase delivery route and a third phase pickup route (Depot-S 3 -S 4 -Depot).

Local search
The solution obtained from the shaking procedure needs to be further improved to obtain a local optimum. The local search procedures are successfully hybrid in the metaheuristics. For example, in Prins (2004) and Vidal et al. (2012), local search procedures were applied to improve the quality of the offspring solutions in the genetic algorithms. Prins (2004) noted that the local search improvement procedure was a key idea of algorithm design. Thus, in our VNS, four kinds of local search methods are used, including interroute and intraroute methods, i.e., interroute 1-0 relocation move, interroute 1-1 exchange move, intraroute 1-0 relocation move and intraroute 1-1 exchange move. Such operators are widely used by VRP studies (Bräysy and Gendreau, 2005;Vidal et al. 2012). The interroute 1-0 move refers to relocating a customer from its current position to a position in another route or constructing a new route, while the intraroute 1-0 move Fig. 2 An example of a saving method Randomly choose a satellite from a first phase delivery route and relocate it into a random position of another first phase delivery route, or construct a new route that only contains this satellite N 2 Randomly choose a satellite from a first phase delivery route and exchange it with one satellite that is randomly selected from another first phase delivery route N 3 Randomly choose two sequential satellites from a first phase delivery route and relocate them into a random position of another first echelon delivery route, or construct a new route that only serves such two satellites N 4 Randomly choose two sequential satellites from a first phase delivery route and exchange them with two satellites that are randomly chosen from another first phase delivery route N 5 Same as N 1 but for the third phase pickup route N 6 Same as N 2 but for the third phase pickup route N 7 Same as N 3 but for the third phase pickup route N 8 Same as N 4 but for the third phase pickup route N 9 Randomly close a satellite. Remove the customers assigned to this satellite one by one; insert each customer into a random position of another second phase delivery and pickup route, or construct a new second phase route that only contains this customer N 10 Randomly choose a second phase delivery and pickup route and remove a random number of sequential customers (bounded by the customer number of the selected route) and insert them to another second phase delivery and pickup route, or construct a new second phase route that only contains such customers C 4 (4,2) C 5 (1,1) C 6 (4,4) C 7 (2,2) first phase delivery route second phase delivery and pickup route Satellite Client Depot S 1 S 2 S 3 C 1 (3,1) C 2 (2,2) C 3 (1,1) C 4 (4,2) C 5 (1,1) C 6 (4,4) C 7 (2,2) Fig. 3 Neighborhood N 1 operator denotes relocating a customer from its current position to another position on the same route. Similarly, the interroute 1-1 move involves exchanging two customers' positions of different routes, and the intraroute 1-1 move exchanges two customers' positions on the same route. The four procedures are used for all three vehicle route phases. Note that for the second phase delivery and pickup route, the interroute 1-0 move may change the status of satellites. For example, now satellite j has only one customer, and satellite k is closed. When this customer is moved from satellite j to satellite k, satellite j becomes closed, and satellite k is open. During the local search procedure, the algorithm adopts the ''first-accept'' strategy, i.e., during a local search procedure once a better solution is found, it is adopted as the new seed for repeating this local search. The local search procedure continues until no improvements can be made.
In terms of the classic local search methods, if a local search starts from a feasible solution seed, the algorithm limits the search in a feasible solution area. Otherwise, if it starts from an infeasible seed, the whole feasible and infeasible solution areas are allowed to be searched. In this paper, the algorithm adopts a new local search strategy. If the solution obtained from shaking (the origin of the whole local search) is feasible, then the local search is restricted to the feasible region, which means that the vehicle capacity cannot be exceeded at any move of the local search. Otherwise, if the shaking result is infeasible, the local search is extended to the infeasible area where the generalized cost, including penalties, is used to evaluate each move.

Acceptance criterion
After the shaking and the local search procedures have been performed, the solution obtained at the current iteration is compared to the incumbent solution to decide whether it is accepted. To prevent the VNS from quickly becoming stuck in a local optimum, inspired by simulated annealing (SA) (Kindervater and Savelsbergh 1997), the scheme of accepting a new but worse solution is adopted. If after shaking and local search procedures, a new solution S 00 that is better than the incumbent solution S is identified, S 00 is accepted directly to replace S. Otherwise, this new solution is accepted with a probability exp Àðcðs 00 Þ À cðsÞÞ ð Þ =T, where c(S 00 ) and c(S) are the generalized costs of two solutions. The annealing temperature T decreases linearly from an initial value T 0 , i.e., after each VNS iteration, T is reduced by T 9 (1/ite), where ite is the VNS maximum iteration time.

Computational results
In this section, computational experiments are conducted to assess the performance of the proposed approach. All algorithms in this paper were coded in C??. The computational experiments were conducted on an Intel E5-2670 processors clocked at 2.6 GHz with 2 GB of memory running Linux. All the algorithms are run 10 times for each test instance. The best and the average solution costs and the average algorithmic running time are obtained from 10 runs for each test instance. C 2 (2,2) C 3 (1,1) C 4 (4,2) C 5 (1,1) C 6 (4,4) C 7 (2,2) S 4 Fig. 4 Neighborhood N 9 operator A variable neighborhood search algorithm with constraint relaxation for the two-echelon… 8889

Test instances from the 2E-VRPSDP benchmarks
Since there is no benchmark instance for our problem, the 2E-VRP benchmark instances from the literature are modified. Thanks to Breunig et al. (2016), the existing 2E-VRP instances were summarized and are available at https://www.univie.ac.at/prolog/research/TwoEVRP. The 2E-VRPSDP test instances are derived from existing 2E-VRP benchmark instances as follows. We select 46 basic 2E-VRP instances in which twenty 2E-VRP instances contain 50 customers, twenty have 100 customers, and the remaining six instances have 200 customers. For each basic 2E-VRP instance, we derive a new instance as follows. We select first 50% customers and set 1/3 of the original demand as the delivery demands and 2/3 of the original demand as the number of pickups. For the left 50% customers, we set 1/3 and 2/3 of the original demand as the pickup and delivery demands, respectively. Other parameters, including locations of all points, vehicle number limits, and vehicle capacity, are the same as those of the original 2E-VRP instance. Since there are only six 2E-VRP instances of 200 customers available in the literature, we change the vehicle capacity and the vehicle number of an echelon and derive another fourteen 2E-VRPSDP instances with 200 customers (with signs of c-g in Table 8). In total, we have sixty 2E-VRPSDP instances for computational experiments. The original 2E-VRP benchmark set, which each of our new 2E-VRPSDP instances comes from, is marked at the beginning of the instance name, and ''SDP'' is added at the end. For example, the new instance generated based on ''E-n51-k5-s2-17'' (2E-VRP benchmark, Set 2b) is named ''Set2b_E-n51-k5-s2-17_SDP.''

Parameter tuning
Based on the above test data, we select 18 2E-VRPSDP test instances to tune the parameters in our algorithm, consisting of six small instances with 50 customers, six medium instances with 100 customers, and six large instances with 200 customers. Each instance is solved by the algorithm with a one-parameter setting ten times. We obtain ten computational results for each test instance, including the best solution cost and mean solution cost among the 10 running results. Then, we calculate the average of all best solution costs and the average of the mean solution cost across all 18 test instances. We use these two statistical results to assess the solution quality of such parameter settings. Table 2 summarizes the final parameter settings of our VNS algorithm used in the experiments. In the following subsections, some parameters are tuned, and a sensitivity analysis is executed to elucidate the effects of the components of the proposed algorithm.

Parameter T 0
Parameter T 0 is the initial temperature of the simulated annealing-based acceptance criterion (Sect. 4.5). Clearly, this parameter has an important effect on the VNS, which controls the probability of accepting a relatively worse solution in the algorithm iterations. We test parameter T 0 on the set of {50, 75, 100, 125, 150, 200} while using the values in Table 2 for other algorithmic parameters. The results are summarized in Table 3. In this table, row ''Avg-Best'' shows the average value of all best solution costs to a total of 18 test instances, and row ''Avg-Mean'' gives the average of all mean solution costs. We find that the solution quality improves with the increase in parameters from a minimum value of 50 and peaks at a value of 100, and then the accuracy of the algorithm decreases with the value of this parameter.

Parameter u
In the VNS, for each solution S, an extended cost c(S) is defined as c(S) = d(S) ? a 9 e(S) for solution evaluation (see Sect. 5.1). At each VNS iteration, when a solution S 00 is obtained after shaking and local search, the algorithm judges its feasibility. If this solution is feasible, the value of a is divided by 1 ? u; otherwise, a is multiplied by 1 ? u. Clearly, this parameter u controls the change frequency of parameter a and the frequency at which the algorithm jumps between feasible and infeasible solution areas. We test u in the set of {1.01, 1.05, 1.1, 1.2, 1.3, 1.5}. The results are shown in Table 4. The best algorithm performance is yielded at u equal to 1.05.

Sensitivity analysis of algorithmic components
Furthermore, a set of experiments is designed to elucidate the effects of VNS algorithm components. We attempt to remove one or a set of components from the algorithm, and the remaining algorithm's performance is explored. Based on these experiments, the roles of the algorithm components are tested and verified. The first is the shaking component (see Sect. 4.3). Ten interroute neighborhood structures (namely, N 1 -N 10 ) are adopted in VNS as shaking. In this part of the experiments, such neighborhood structures are classified into four sets. The first set contains N 1 , N 3 , N 5 , and N 7 , which relocate a satellite to another position; the second set consists of N 2 , N 4 , N 6 , and N 8 , which exchange the positions of two satellites. The third is N 9 , which closes a satellite, and the final fourth is N 10 , which removes a number of sequential customers and inserts them into another route. We sequentially delete one set of neighborhood structures from the whole algorithm and yield four versions of algorithms, i.e., ''no-relocation,'' ''no-exchange,'' ''no-close,'' and ''no-remove.'' We also aim to test the component of the local search (Sect. 4.4). The solution obtained from the shaking procedure is further improved by a set of local search methods. The first set of local search methods is 1-0 relocation, and the second is 1-1 exchange (details are referenced in Sect. 4.4). Again, each time one set of local search methods is shielded, the remaining algorithms ''no-1-0 LS'' and ''no-1-1 LS'' are tested. Each version of the modified algorithm is executed 10 times over each instance. The results are shown in Table 5. The results show that all these algorithmic components enhance the performance of the algorithm because ''Whole VNS'' yields the best results. Therefore, it is concluded that the algorithm performance will be compromised if one component is deleted from the algorithm.
Furthermore, in terms of each of the above operators, we also want to know how many times it is used, or whether it truly improves the solution during the VNS iteration process. For example, considering four types of Shaking operators (relocation, exchange, close, and remove), we record the mean execution times of each operator during VNS iterations. They contribute 12.7%, 10.7%, 37.2%, and 39.4% of shaking executions, respectively. This statistical result shows that each shaking operator is executed and has a considerable positive effect on the algorithm performance. We also record the success rate of two local search methods 1-0 relocation and 1-1 exchange (i.e., the number of times that one LS operator improves the solution during VNS whole iterations, over the total improvement times of all LS operations). It is found that 1-0 relocation contributes approximately 53% local search improvement, whereas 1-1 exchange contributes 47% improvement. This result also indicates that two local search methods are necessary in the algorithm. Furthermore, we find that the most running time of VNS is used by local search procedures, and they account for approximately 90% of the running time of the whole algorithm. The remaining time is mainly used by shaking operators. Furthermore, the computational time of each shaking operator accounts for a similar percentage of the total shaking time as their executions (i.e., approximately 10%, 10%, 35%, and 45%). Similarly, each local search operator's running times also account for a similar percentage of the total consumption time as their contributions stated above.

Computational results for 2E-VRPSDP instances
We cannot find other algorithms for the 2E-VRPSDP, so in this paper, we first develop a simplified version of our approach, named VNS1, to test the design and performance of the proposed VNS. In the VNS1 approach, the basic framework of the algorithm is the same as that of the above approach, including the construction of an initial solution, shaking, local search and the acceptance criterion. However, the search space of VNS1 is always restricted to the feasible region, which means that in the shaking and local search procedure, the vehicle capacity constraints are always obeyed. Thus, the solution obtained at any iteration  of VNS1 is always feasible, and its penalty cost equals zero. Such a version of VNS is more classic and standard. The VNS approach proposed above (extending the algorithm to an infeasible solution area) is denoted as VNS2.
To give a fair comparison between two approaches, VNS1 and VNS2, for each test instance, we change the stop criterion of algorithm VNS1, i.e., we do not use maximum iterations to terminate VNS1; it is stopped after a special total running time that equals the computational time of VNS2. That is, the two approaches have the same running times for solving each test instance. Additionally, we also try to solve the problem using the commercial solver Gurobi and a decomposition-based approach. In terms of the Gurobi solver, we implement the mathematical model in Sect. 3 and solve each test instance with a maximum of 3 h running time much longer than that of the VNS. The decomposition-based approach is built on the basis that the 2E-VRPSDP is decomposed into two levels, in each of which a special vehicle routing problem is solved. That is to say, the classical VRPs are solved at the first level, and a pickup and delivery problem is optimized at the second level (refer to Fig. 1). Finally, the solutions of two levels are combined to obtain a whole solution to our 2E-VRPSDP problem. For the routing problem at each level, a simplified VNS algorithm of this paper is used to optimize the vehicle routes. The details of the simplified VNS algorithms are presented in Appendix.
We tested a total of sixty 2E-VRPSDP instances on Gurobi, decomposition-based approach, VNS1 and VNS2. The numerical results are aggregated in Table 6, according to the scales of the test instances. The detailed results are presented in Table A1-A3 in Appendix. In Table 6, columns ''Avg'' and ''Best'' give the average solution cost and best solution cost over ten runs of the corresponding algorithm, respectively.
As shown in Table 6 and Tables A1-A3 in Appendix, for all instances, VNS2 yields better solutions from the perspective of solution cost (quality). VNS2 can find a better solution for each test instance than Gurobi, decomposition-based approach, and VNS1. As shown in Table 6, compared with Gurobi, VNS2 can improve mean solution costs quality 17.72%, 27.77%, and 40.43% for 50-customer, 100-customer, and 200-customer instances, respectively. The superiority of our VNS2 over the solver becomes clearer for larger-scale tests. Compared with the decomposition-based approach, the superiority of our VNS is still clear: For three scales of test instance, the mean solution costs of our VNS2 is smaller than that of the decomposition-based approach 25.38%, 15.00%, and 12.42%.

Computational results for 2E-VRP benchmarks
To further assess the proposed VNS approach, we tested it on existing 2E-VRP benchmark instances from the literature. The 2E-VRP can be regarded as a special 2E-VRPSDP when the pickup demand of each customer of the 2E-VRPSDP is zero. In such cases, the second phase delivery and pickup route is reduced to the second phase delivery route; the third phase pickup routes do not exist. Note that the classic 2E-VRP assumes that each satellite can be visited by more than one first phase delivery route; however, in the 2E-VRPSDP, a satellite is limited to service from no more than one vehicle. Thus, we select twenty  (2014) is carried out using a single core of an Intel Pentium Dual-Core E5500 processor 2.8 GHz and 2 GB of memory. The PLNS in Mühlbauer and Fontaine (2021) is performed on a computer with 8G RAM and an Intel i5-6200 processor with 2.4 GHz, two cores and four threads. The PLNS is a parallelized algorithm, which is executed on the four threads in parallel. In addition to the results of the above algorithms, we also compare our results with the current BKS of each test instance from the literature. Table 7 aggregates and compares detailed computational results on the existing 2E-VRP benchmark instances, according to the scale of test instance (customer number). The detailed results for each test instance are shown in Table A4 in Appendix. In Table 7, Column ''BKS'' refers to the best-known solution of the instance from the existing literature. For the other approaches (LNS-2E, ALNS, GRASP ? VND, PLNS), subcolumn ''Best'' gives the best objective value found within their five runs; column ''Avg'' shows the average objective value of these five runs. Note that Zeng et al. (2014) only provided the average running time (in seconds) to find the best solutions (not the whole algorithm running time) in their paper. Therefore, for their results, Table 7 gives such times, whereas for other approaches, the whole algorithm running times (average run one time, in seconds) for solving each instance are shown.
First, we compare the solution quality of different algorithms. As shown in Table 7 and Table A4, the first four approaches can find BKS for the first ten instances that contain 50 customers. The VNS algorithm and the LNS-2E of Breunig et al. (2016) have better performances since they can obtain BKS for every algorithm running, i.e., the average solution (column Avg) equals the best solution (column Best). Algorithm PLNS in Mühlbauer and Fontaine (2021) fails to find the BKS for three test instances (Instance50-2, Instance50-37, and Instance50-49).
Considering larger-scale test instances (100 customers), we observe that VNS has the highest quality mean of the best solution (1128.42). Considering the largest instances (200 customers), the mean of the VNS best solution (1372.73) is only worse than that of the PLNS (1355.27) and is better than the other algorithms. Regarding the number of BKS found by each approach, as shown in Table A4, VNS finds five existing BKS; LNS-2E and PLNS both obtain three existing BKS. Compared with LNS-2E, VNS can find six better solutions, and LNS-2E can only find one better solution than VNS. Compared with PLNS, VNS finds three better solutions (100-5-1, 100-5-1b, and B-n101-5), and for 3 test instances, PLNS obtains better solutions (200-10-1, 200-10-1b, and 200-10-2). These results indicate that for such large test instances, VNS is better than LNS-2E, and the two algorithms VNS and PLNS, are similar in solution quality.
Additionally, Table 8 aggregates the best objective value found by each approach during all experiments, including those used for parameter calibration. We can observe the superiority of the VNS and PLNS over other approaches. The detailed results of Table 8 are shown in Table A5 in Appendix. As shown in Table A5, for 15 out of a total of 20 test instances, the VNS can find the BKS from the literature. The BKS found by each other algorithm is much less than VNS. Such numerical results indicate that the VNS can obtain high-quality solutions for the twoechelon vehicle routing problem with simultaneous delivery and pickup.
In addition to the solution quality, in terms of runtime for such test instances, as shown in Table 7, the mean running time of VNS is shorter than that of LNS-2E and ALNS. For example, the mean computational times of VNS and LNS-2E over all test instances are 284 s and 480 s, respectively. Note that the total runtime of GRASP ? VND is not given in the literature; its running time shown in Table 7 is only the running time to find the best solutions (not the whole algorithm running time). We also found that the running time of PLNS is shorter than that of our VNS, especially for large-scale test instances with 100 ? customers. However, please note that the PLNS is a parallelized metaheuristic, and our VNS is not designed as a parallel algorithm. The parallelized technique can greatly reduce the running time of PLNS. In summary, compared with existing state-of-the-art approaches, the VNS in this paper does not have an obvious superiority in computation time. However, its running time is acceptable and reasonable for all these test instances.

Conclusions and future work
In this paper, a new two-echelon transportation problem was investigated, the two-echelon vehicle routing problem with simultaneous deliveries and pickups (2E-VRPSDP). The 2E-VRPSDP considers the distribution and collection of freight on two echelons. It is more complex than the classic two-echelon vehicle routing problem (2E-VRP) and the vehicle routing problem with simultaneous deliveries and pickups (VRPSDP). To address this challenging problem, a VNS algorithm was developed, which extends the search space into the infeasible region and strikes a balance between infeasible solutions and feasible solutions by adding a linear penalty cost function. To prevent the algorithm from becoming stuck in a local optimum, an acceptance criterion inspired by simulated annealing was embedded into the proposed method. To examine the algorithmic performances, new 2E-VRPSDP test instances were generated based on the existing 2E-VRP benchmarks. Our approach was compared with its simplified version, which restricts the search in the feasible region. The results show the priority of the proposed method in obtaining higher-quality solutions. The proposed approach was also compared with existing state-of-the-art algorithms on the 2E-VRP benchmarks from the literature. The numerical results prove the efficiency and effectiveness of the proposed approach.
The work presented in this paper can be extended in two directions. The first is that uncertain elements can be added to the current framework, such as customers with random delivery and pickup demands. Stochastic demands may make the problem more similar to real-life applications. Additionally, this paper focuses on a heuristic algorithm, and the optimal solution cannot be guaranteed. It is interesting to consider exact decomposition algorithms such as bender decomposition or branch-and-price algorithms to solve this problem optimally in future work.