An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation

ARTICLE IN PRESS

JID: EOR

[m5G;May 2, 2017;13:24]

European Journal of Operational Research 0 0 0 (2017) 1–15

Contents lists available at ScienceDirect

European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor

Production, Manufacturing and Logistics

An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation Farshid Azadian a,∗, Alper Murat b, Ratna Babu Chinnam b a b

College of Business, Embry-Riddle Aeronautical University, 600 S. Clyde Morris Blvd., Daytona Beach, FL 32114, USA Department of Industrial and Systems Engineering, Wayne State University, 4815 Fourth Street, Detroit, MI, 48202, USA

a r t i c l e

i n f o

Article history: Received 11 July 2016 Accepted 16 April 2017 Available online xxx Keywords: Transportation Pickup and delivery problems Air cargo Multi-airport region

a b s t r a c t This study considers a freight forwarder’s operational implementation of alternative access airport policy in a multi-airport region for cargo transportation. Given a set of heterogeneous air cargo customers, the forwarder’s problem is to simultaneously select air cargo flight itineraries and schedule the pickup and delivery of customer loads to the airport(s). This problem is formulated as a novel pickup and delivery problem, where the delivery cost is both destination and time dependent. A mixed-integer linear model in presented and an efficient solution method based on decomposition is developed. The performance of the solution algorithm is evaluated by computational experimental study. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Over the past two decades, air mode of transportation for carrying goods around the globe has become a vital part of many companies’ logistics. The demand for air cargo demonstrated steady growth and is expected to maintain the trend despite the economic crisis of the past few years. Annual forecast reports by Airbus (2015) and Boeing (2015) predict 4.4 to 4.7 percent annual growth for global air cargo tonnage over the next 20 years, respectively. Similar to a growing demand for air transportation, although not at the same rate, the aviation industry has also expanded and air transportation has become more accessible to shippers with more options and at more locations. In the domestic US market, construction of new airports, capacity expansion of existing airports, conversion of some military airports to civilian airports, increase in the aviation network connectivity, availability of more flight itineraries, and improved road accessibility of airports cause the service zones of airports to expand and, often, overlap. This has resulted in the creation of Multi-Airport Regions (MARs) where several airports are accessible in a region and can substitute and supplement each other in meeting the region’s demand for air transportation (Loo, 2008). Hall (2002) proposed the Alternative Access Airport Policy (AAAP) and argued that for air cargo transportation considering multiple airports (and subsequently flight itinerary options) in a MAR can be beneficial

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (F. Azadian), [email protected] (A. Murat), [email protected] (R.B. Chinnam).

to reduce truck mileage, decrease sorting and handling costs, improve delivery service level, and avoid congestion on both road and air networks. He has used the case of Southern California to point out the potentials for AAAP in addressing the complexities of responding to increasing demand in that region. While Hall (2002) outlined and discussed the advantages of the AAAP, to the best of our knowledge, there is no study on the modeling and implementation of AAAP. In this paper, we consider a freight forwarder’s operational implementation of AAAP in air cargo transportation. The freight forwarders, who are intermediaries between shippers and carriers, are responsible for more than 90% of air cargo shipments (Hellermann, 2006). Shippers are often only concerned with the safe arrival of their cargo to the destination at a reasonable price and increasingly move toward leaving all the routing and booking decisions to freight forwarders. Accordingly, freight forwarders are responsible for finding the best flight itinerary options and booking capacity with air carriers to get the cargo to its destination before delivery due date. A flight itinerary option refers to a set of one or multiple connecting flights that carry cargo from its origin to its final destination. Various flight itinerary options may share some flight legs. For instance, a flight leg from the origin airport to an airline hub is shared among all the flight itinerary options that fly through that hub. The process of routing cargo on a scheduled air network to establish the best flight itinerary options and calculate their costs and destination arrival times is a complex problem. Azadian, Murat, and Chinnam (2012) studied this problem and provided an algorithm for air cargo routing and dynamic rerouting. We rely on their research for identifying the flight itinerary options and calculating their costs. In the case of a MAR and under the AAAP, there

http://dx.doi.org/10.1016/j.ejor.2017.04.033 0377-2217/© 2017 Elsevier B.V. All rights reserved.

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

JID: EOR 2

ARTICLE IN PRESS

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

may be a large number of feasible flight itineraries that need to be considered for each customer’s load (see Appendix C for an illustrative example). In addition to flight booking, forwarders are also responsible for the pickup and delivery of customers’ loads, which is often performed by a fleet of hired trucks. When customers’ loads are lessthan-truck-load, which is a common case for air cargo, consolidating loads of multiple customers in a single truck trip provides an opportunity to reduce logistics costs. In fact, the profitability of a forwarder relies directly on its overall cost efficiency, which depends on the flight itinerary decisions as well as pickup and delivery operations. The combination of these decisions constitutes a very complex operational management problem. Our objective, in this paper, is to introduce and model this practical problem and provide a novel solution methodology for it from the perspective of a freight forwarder. The forwarders problem constitutes of two sets of decisions. One decision component in this problem is the flight itinerary assignment to different customers’ loads that are geographically dispersed in the MAR. These decisions are driven by the availability of flight itinerary options, cargo acceptance cut-off times, destination arrival times, flight itinerary costs, and tardiness penalties. The other decision component is the multi-vehicle routing to pick up customer loads and deliver them to the airports prior to the cut-off times for cargo acceptance of the assigned flight itineraries. These routing decisions are affected by the locations (depot, customers, and airports), starting times of the selected flight itineraries, and the vehicle fleet size. The structure of the forwarders problems generalizes the Pickup and Delivery Problems (PDP) in several aspects and demand the introduction of a new type of the PDPs. Our contributions in this research are threefold. First, we introduce and model the operational implementation of the AAAP for freight forwarders in a MAR, which generalizes several known pickup and delivery problems in terms of model structure and objective function. Second, we propose a novel and efficient decomposition based solution methodology, which can also be used for the new class of pickup and delivery problems. Finally, we present the results of an experimental study to evaluate the performance of the proposed algorithm. The rest of this paper is organized as follows. We briefly review the relevant literature in Section 2 and point out the gap in the existing literature. In Section 3, we present a detailed account of the problem modeling and formulation. In Section 4, we demonstrate the shortcomings of traditional methods in solving the problem and develop a solution methodology that is able to solve the problem efficiently. In Section 5, we evaluate the performance of our solution method and report on the results of the computational study with experimental problem instances. Section 6 concludes the paper with a discussion and future research directions. 2. Related Literature The freight forwarder’s operational implementation of the AAAP is closely related to the pickup and delivery problem. Pickup and delivery problems have been extensively studied in the past decades; for a comprehensive survey see Berbeglia, Cordeau, and Laporte (2010), Berbeglia, Cordeau, Gribkovskaia, and Laporte (2007), Laporte (1992, 2009), Parragh, Doerner, and Hartl (2008a, b), and Toth and Vigo (2001). Generally, the PDP involves routing a fleet of vehicles to satisfy a set of transportation requests between the given origins and destinations. In the PDP, all the origin pickups must precede the destination deliveries and be performed by the same vehicle. Moreover, each route must start and terminate at the same location (i.e., depot). The PDP usually considers capacitated vehicles and the goal is to minimize a criterion related to a travel measure. The travel measure can be as simple as the total

travel distance for urban commercial vehicles (Figliozzi, 2007) or more complex as the total excess riding time over the direct ride time in passenger transportation (Diana & Dessouky, 2004). The PDP can be classified into two categories: transportation between customers and the depot, and transportation between the pickup and delivery locations (Parragh et al., 2008a). The proposed problem in this paper falls into the latter category, as the depot (the starting point of trucks) does not generate or accept any load. This category can be further classified into paired and unpaired pickup and delivery locations. In the paired PDP, also known as the One-to-One PDP (1-1-PDP), the load picked up from a customer location can only be delivered to one of the delivery locations associated with the given customer. Some customers, however, may share the same delivery location. In the stacker-crane problem (SCP), unit loads of nonidentical commodities have to be transported from the origin to the destination using a unit capacity vehicle (Frederickson, 1978). In the Vehicle Routing Problem with Pickup and Delivery (VRPPD), the unit capacity requirement of SCP is relaxed and replaced with a set of constraints based on the load properties (e.g., weight, volume, or unit count). A special case of the VRPPD is the VRPPD with Time Windows (VRPPDTW) where visiting a pickup or delivery location is only allowed during a predefined time window. While the VRPPD generally concerns goods transportation, the dial-a-ride problem (DARP) addresses passenger transportation and therefore includes additional side constraints (e.g., maximum ride time, time windows, or service quality). Moreover, the objective function optimizes measures such as customers (in)convenience; see Cordeau and Laporte (2007), Kirchler and Calvo (2013), and Paquette, Cordeau, Laporte, and Pascoal (2013) for a survey of recent modeling and solution algorithms for DARP. In comparison, the unpaired PDPs, also known as Many-to-Many PDP problems (M-M-PDP), consider the case where any commodity can be picked up and delivered to delivery locations that accept the commodity. The M-M-PDP was initiated with Anily and Hassin (1992) that introduced the swapping problem (SP) for moving ncommodity objects between customers with a single unit capacity vehicle. In the SP, each customer supplies one type of commodity and demands a different type. In addition to the n-commodity case of the SP, there are several other single commodity problems that are studied under the M-M-PDP where picked up loads are homogenous. Hernández-Pérez and Salazar-González (2004a, b, 2007) introduced and studied the one-commodity pickup and delivery traveling salesman problem (1-PDTSP). The 1-PDTSP is a more general case of the Q-delivery traveling salesman problem (Q-DTSP) by Chalasani and Motwani (1999) and the capacitated traveling salesman problem with pickup and deliveries (CTSPPD) by Anily and Bramel (1999). In the 1-PDTSP, a single vehicle, starting from a depot, transports goods from pickup nodes to delivery nodes without exceeding the vehicle capacity; the objective is to minimize the total traveling cost. Q-DTSP and CTSPPD are special cases of 1-PDTSP where the pickup and delivery quantities are all one unit and the vehicle capacity is restricted (i.e., Q units). Hernández-Pérez and Salazar-González (2009) later extend their 1-PDTSP to the MultiCommodity One-to-One Pickup and Delivery Traveling Salesman Problem (m-PDTSP); however, with this extension, the problem is not an M-M-PDP anymore. The multi-commodity variation was studied by Hernández-Pérez and Salazar-González (2014). In the multicommodity variation, a given picked up commodity can be used to satisfy any demand of that commodity disregarding the origin. The problem of the forwarder’s operational implementation of the AAAP is essentially a PDP as it consists of transporting loads from customer sites (pickup locations) to the airports (delivery locations) in a MAR. The depot is both the origin and final destination of the vehicles, although, it is neither pickup nor a delivery point. Despite the similarities, the proposed problem has

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

several major differences with existing PDPs and it cannot be cast as equivalent to any one of the problem classes. The major source of distinction of this problem is the assignment decision of customers to flights, which in turn impose airport delivery time constraints and affect the objective function through the inclusion of flight assignment costs The lack of a predefined assignment between customers and flight itineraries is what separates the proposed problem from 11 PDPs. In other words, if flight itinerary assignments were fixed (e.g., there were only one flight option to assign for each customer out of all the airports), then this problem (with some transformation and tricks) can be presented as 1-1 PDP with time windows (i.e., VRPPDTW). Indeed, under this scenario, even the problem of confirming that a feasible solution exists is NP-complete (Savelsbergh, 1985). Of course, the implementation of AAAP consists of a large number of options for flight itinerary assignments, resulting in further combinatorial complexity. The proposed problem differs from M-M PDPs in that the delivery cost of customer loads is both time and destination dependent. Similar to the M-M PDP, in this problem, a customer load can generally be delivered to any of the airports in a MAR. However, considering the flight itineraries offered by different airports, some of the airports might never be visited. Furthermore, delivering a customer’s load to different airports results in different costs based on the available flight itinerary options at the time of delivery. As can be seen, aside from the time windows, the assignment decision of customers’ loads to flight itinerary options distinguishes this problem from the existing PDPs, in which the objective is limited to optimizing the total travel distance/time or a measure related to it. Due to this extension to the PDP, this problem cannot be transformed to any existing PDP and it requires a new solution approach. We denote this problem as the PDP with Assignment and TimeDependent delivery cost (ATD-PDP). The use of the term assignment indicates that the delivery cost of a customer’s load depends on the airport and flight itinerary selected. In our problem, time dependency refers to the fact that flight itinerary options at each airport have a given cargo acceptance cut-off time; consequently, the flight itinerary options (and in turn the cost of transporting customer’s load) changes by time at each airport. The proposed problem characteristics have not been studied in the literature and, to the best of our knowledge, this is the first research on PDPs with assignment and time dependent delivery costs. The proposed problem is clearly an NP-hard problem in the strong sense as it coincides with the VRPPD when there are only one airport and one single itinerary (accepted by all customers), which departs late enough to complete all pickups and deliveries. 3. Problem description and formulation In this section, we develop the model formulation of the ATDPDP. We first briefly define our problem and discuss the time dependent delivery cost, which is a critical component of the objective function as it jointly captures the flight itinerary costs and tardiness penalties. Next, we describe the graph transformation and present the mixed integer programming model formulation. Last, we introduce and discuss pre-processing steps and valid inequalities to strengthen the formulation. Our problem concerns a freight forwarder’s operational planning of the pickup and delivery of different customers’ air cargo loads that are geographically dispersed via a fleet of trucks and delivering to different airports in a region based on the flight itinerary options to be selected for each customer air cargo load. The two sets of decisions are the flight itinerary selection decisions and pickup and delivery decisions. The constraints are time windows due to the departure cut-off time for the flights. The

3

objective is to minimize the total cost consisting of flight itinerary costs, tardiness penalties subject to the destination arrival times (i.e., time dependent delivery cost), and road transportation costs. We assume fixed number of identical vehicles without any capacity constraint. We further assume no splitting of the pickup and delivery as it is not in practiced in the air cargo transportation. Let Go = (Vo, Eo ) be a directed graph representing the network topology of the problem where Vo is the set of nodes and Eo is the set of connecting edges. The set Vo consists of the depot d, the set of customers (pickup locations) C, and the set of airports (delivery locations) H; i.e., Vo = {d} ∪ C ∪ H. Let K be the set of uncapacitated homogeneous vehicles (trucks) that originate from the depot and op operate during the depot’s opening (θd ) and closing hours (θdcl ). A cost ci j and a travel time ti j is associated with each edge ∀(i, j ) ∈ Eo, i = j of the network, where ci j ≥ 0 and ti j ≥ 0. We assume that the triangle property holds for the travel times and travel costs; i.e., we have ti j + t jg ≥ tig and ci j + c jg ≥ cig , ∀i, j, g ∈ Vo. Note that, if needed, the fixed cost of utilizing a vehicle can be captured by adjusting the cost parameter cdj , ∀ j ∈ C ∪ H. We further assume that travel time ti j is deterministic and time-independent. We assume that there are no time windows for customers’ pickups and the service times (e.g., loading and unloading) are negligible. The formulation can be easily extended to incorporate these considerations, as the methods presented do not rely on their absence. Moreover, due to the nature of air cargo we do not allow transshipment and/or split pickup or delivery. Let Rh be the set of flight itinerary options available at airport h ∈ H on the day of operation. The flight itinerary options consist of the outgoing flights from the given airport and subsequent connecting flights that reach to customers destinations. This set can be pre-identified for each individual customer by evaluation of the available options considering the customer load requirements, destinations and time-dependent delivery costs (i.e., prohibiting excessive transportation duration). Considering the connected air network, there is often more than one flight itinerary option for each customer to fulfill the shipping request. The cost of assigning a flight itinerary r ∈ Rh to customer i is Firh , which accounts for the transportation cost of the carrier as well as the delivery service level related costs, such as tardiness penalties. The starting time of the flight itinerary r is denoted by Qrh (i.e., the cargo drop-off cut-off time for the first flight of the itinerary). We only consider those flights that can be op used on the day of operation, e.g., Qrh ≥ θd . As mentioned before, a flight itinerary option is a direct flight or combination of connecting flights that can be used to carry a customer load to its destination. Indeed, different flight itinerary options that end in different destinations may share the same first leg. Consequently, those given flight itinerary options have the same cut-off time (Qrh ) but result in a different cost (Firh ) for each customer, considering the destination and individual customer load characteristics. If a specific customer load cannot be assigned to a given flight itinerary option available at the regional airports, that option can be eliminated from further consideration for that given customer. We will comment on this scenario in Section 3.3. 3.1. Time dependent delivery cost In assigning the customer i’s cargo to a flight itinerary r ∈ Rh , the freight forwarder accounts for the airport h arrival time. The assignment is feasible only if the airport delivery time (t) is on or before the flight itinerary starting time, i.e., t ≤ Qrh . When a customer’s load is delivered to an airport h at time t and there are no flights available, t > maxr∈Rh {Qrh }, then the air cargo is assigned to a recourse flight itinerary r0 ∈ / Rh , e.g., a next day itinerary. We assign a penalty cost Fih0 > Firh ∀r ∈ {Rh , h ∈ H }, for airport delivery after the departure time of the last flight on the day of operation. Accordingly, we define the time dependent airport delivery cost of

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR 4

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

Fig. 1. Illustrative airport h delivery cost function for customers i, j ∈ C; customer i has two flight itinerary options (left) and customer j has a single flight itinerary option (right).

delivering customer i’s load to airport h at time t, f (h, i, t ), as follows:



f (h, i, t ) =

min{Firh |t ≤ Qrh } r∈Rh

Fi0h ,



i f t ≤ max Qrh



r∈Rh

otherwise

As implied in this definition, we do not need to consider all itinerary options for each customer. We can identify the potential set of itinerary options that are dominated by at least another itinerary option from the same airport. The flight itineraries that are dominated for all customers are removed from further consideration. The flight itineraries that are dominated only for a subset of customers are preprocessed such that their assignment to that subset of customers is precluded. Lemma 1 provides the conditions necessary to identify the dominated flight itineraries from airport h for customer i. Lemma 1. Given two flight itineraries r, r ∈ Rh , r = r , itinerary r is dominated by itinerary r if (a) Firh ≤ Firh and Qrh < Qrh or (b) Firh <

Firh and Qrh ≤ Qrh . Moreover, if (c) Firh = Firh and Qrh = Qrh , considering either one is sufficient. Proof. The proof is evident from the definition of f (h, i, t ).



Upon the elimination of dominated itineraries, the following corollary states that there exist no two flight itineraries for customer i at airport h that either depart at the same time or have the same cost. Corollary 1. After eliminating the dominated flight itineraries, there are no two flight itineraries such that r, r ∈ Rh , r = r in f (h, i, t ) h h h h with Qr = Qr or Fir = Fir . Theorem 1 characterizes the airport delivery cost function after eliminating the dominated itineraries. Theorem 1. The airport delivery cost function f (h, i, t ) based on non-dominated flight itineraries is a non-decreasing step function with discontinuities at every Qrh ∀r ∈ Rh . Proof. Let us first consider the single flight itinerary case where f (h, i, t ) = Firh , ∀t ≤ Qrh and Fih0 otherwise. Since Fih0 > Firh , the cost function f (h, i, t ) is a step-function, which is non-decreasing and has a single discontinuity at Qrh . In the case of more than one flight itinerary, let us consider any two itineraries r, r ∈ Rh . From Lemma 1 and Corollary 1, if Qrh < Qrh then Firh < Firh . Therefore, for any two delivery times t1 and t2 where t1 < t2 , we have f (h, i, t1 ) ≤ f (h, i, t2 ). In this case, f (h, i, t ) is a nondecreasing step-function with discontinuities at Qrh and Qrh . The case for more than two itineraries follows from the induction.

Thus, the airport delivery cost function is a non-decreasing stepfunction with discontinuities at the starting times of the nondominated itineraries.  Fig. 1 illustrates a typical airport delivery cost function at airport h for two customers i, j ∈ C. There are two flight itinerary options available r = 1 and 2. While customer i can use both r = 1 and 2, customer j can only use the flight itinerary r = 1 and its load cannot be shipped by itinerary r = 2, e.g., destination of itinerary r = 2 is different from the customer j’s destination. Note that airport delivery after Q2h for customer i (Q1h for customer j) h for customer j). will result in the penalty cost of Fih0 (Fj0 From Theorem 1, we can deduce that waiting at any node or delaying any airport delivery is suboptimal for ATD-PDP. Moreover, all used vehicles start at their earliest time from the depot. 3.2. Graph Transformation In ATD-PDP, each airport can be visited multiple times by a vehicle to deliver loads from different customers. Consequently, a feasible solution may not be a Hamiltonian cycle. Hence, in modeling the ATD-PDP, we need to keep track of the order of these airport visits for each vehicle by introducing additional variables. Moreover, another set of additional variables is needed to handle the step-function characteristic of the airport delivery cost. To reduce the complexity of the ATD-PDP’s formulation and eliminate the need for these additional variables, we perform a graph transformation of the original network graph Go = (Vo, Eo ). We now describe important properties of the optimal solutions of ATD-PDP used in the graph transformation. The first property relates to preemption, which is the act of temporarily leaving the previously picked up load at a location that is not its destination for retrieving it for delivery at a later time. Lemma 2. There is an optimal solution for ATD-PDP that is not dominated by a solution with preemption. Proof. First, vehicles have no capacity restrictions to motivate preemptive solutions. In addition, a preemptive solution potentially prolongs deliveries by introducing additional node visits. For any solution with preemption, we can identify a similar solution without preemption where the return visit for picking up the dropped load is eliminated while the remainder of the decisions remains the same. Since this elimination does not increase the airport arrival time, then, from Theorem 1, the non-preemptive solution has same or better objective function than that of the preemptive solution. 

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

JID: EOR

ARTICLE IN PRESS

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

5

Fig. 2. Illustration of a sample feasible solution in the original (a) and transformed (b) graphs.

We can deduce for Lemma 1 that in ATD-PDP, there is an optimal solution where all customer nodes are visited at most once. Thus, we can restrict the visit of each customer to at most once. The airport nodes, in contrast, can be visited more than once by each vehicle. However, the following theorem establishes that each vehicle visits an airport only once for each itinerary. Theorem 2. There exists an optimal solution of ATD-PDP where each vehicle delivers customers’ load to an airport for each flight itinerary only once. Proof. Consider an optimal solution in which a set of customers (S) are assigned to a given flight itinerary r at airport h. Assume that these customers are delivered to the airport by one vehicle but in two visits at times t1 and t2 consecutively, where t1 < t2 . Clearly, t1 < t2 ≤ Qrh . Let us denote the set of delivered customers at each visit as two distinctive and non-empty sets of S1 and S2 respectively. To prove the theorem it is sufficient to show that moving all the customers in set S1 to set S2 will result in a feasible solution with the objective value the same as the optimal objective value. First, since set S2 is not empty and vehicles are uncapacitated, the proposed solution is feasible. Moreover, since the same itinerary is used, the objective value is the same as the original optimal value. In other words, although the vehicle may still visit the airport at time t1 for other itineraries, since S1 is empty in the proposed solution, itinerary r is used only once and in the second visit.  Theorem 2 states that we can restrict the solution of ATDPDP to those solutions where each flight itinerary requires at most one visit to the airport. In other words, we only need to consider visits to an airport h equal to the number of flight itineraries ∀r ∈ Rh plus an additional visit for the recourse flight r0 ∈ / Rh . We use this property to perform the graph transformation. In our graph transformation scheme, we partition each airport node h into |Rh | + 1 nodes, each node representing a single flight itinerary. In the remainder, we refer to these nodes as flight nodes. Let G = (V, E ) be the transformed graph of the original graph Go = (Vo, Eo ). In this transformation, each airport node h ∈ H is replaced by |Rh | + 1 flight nodes: |Rh | nodes each corresponding to a flight itinerary plus another node for the recourse flight. Consequently, the airport set H is replaced with a new set of  flight nodes r ∈ R, where |R| = h∈H |Rh | + |H |. The geographical locations of the flight nodes are identical to that of their respective airport nodes. Then, we have V = {d} ∪ C ∪ R. The cost of assigning flight itinerary r ∈ Rh to customer i (Firh ) is replaced with the delivery cost (Fir ) to flight node r ∈ R. Note that we are using the

same index r for itineraries and flight nodes. Further, we introduce a hard upper time window Qr for flight node r, i.e., it cannot be visited after Qr . The flight node for recourse flights has the delivery cost of Fi0 and upper time window of infinity. As for the edges, we replace the airport ∀h ∈ H edges ∀( j, h ) ∈ Eo, ∀ j ∈ Vo\{h} with new flight node edges ( j, r ) ∈, ∀ j ∈ V, ∀r ∈ Rh and assign edge travel times t jr = t jh and costs c jr = c jh . A similar procedure is repeated for the outgoing links. In addition, a new set of links interconnecting the flight nodes are added with zero travel time and cost for the flight nodes generated from the same airport. The transformed graph G = (V, E ) inherits all the edges connecting depot to customers and customers to customers. While a feasible solution in the original graph may not be a Hamiltonian cycle, the same solution is represented with one or more Hamiltonian cycles on the sub-graphs of the transformed graph. Indeed, any solution in graph G can be easily transferred back to a solution in original graph Go by collapsing the flight nodes back to their original airport node. Fig. 2 illustrates the graph transformation on a network with 5 customers and 2 airports, each with 2 flights. In the feasible solution illustrated in Fig. 2a, loads from customer(s) {1},{2, 3, 4}, and {5} are assigned to flight itineraries r2 at airport H1, r3 at airport H2, and r4 at airport H2, respectively. While vehicle 1 s trip is a Hamiltonian cycle, vehicle 2 visits the airport H2 twice. In the transformed graph in Fig. 2b, this solution is represented in a single Hamiltonian cycle where vehicle 1 visits flight node r1, and vehicle 2 first visits flight node r3 and then subsequently r4. In Fig. 2b, the shaded flight nodes correspond to flight nodes for recourse flights. The graph transformation eliminates the need for introducing additional integer variables for tracking the order of vehicle visits to airports as well as handling the step-function characteristic of the time dependent delivery cost. This transformation further reduces the complexity of the ATD-PDP’s formulation. In particular, it allows network preprocessing and introducing valid inequalities to strengthen the formulation as described in Section 3.4. 3.3. Problem Formulation The objective of the ATD-PDP is to pick up all customer loads, assign loads to flight itineraries and deliver loads to the airports on time while minimizing the total cost. We now formulate the ATD-PDP using the transformed graph as a mixed-integer programming model. Let xki j denote the binary decision variable indicating whether vehicle k travels from node i directly to node j. Let ykir be the binary decision variable indicating whether the load

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR 6

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15 Table 1 Notation used in our problem formulation. C R d V K

set of customers; indexed by i and j set of flight itineraries; indexed by r sepot, the starting and ending node for all vehicles set of vertices; that is, V = {d} ∪ C ∪ R set of vehicles; indexed by k the depot opening and closing time, respectively acceptance cut-off time of flight itinerary option r ∈ R cost of flight itinerary option r ∈ R for customer i ∈ C travel time from node i to j, i, j ∈ V cost of traveling from node i to j, i, j ∈ V decision variable indicating whether vehicle k ∈ K went directly from node i to j, i, j ∈ V decision variable indicating whether vehicle k ∈ K picked up customer i ∈ C load to be delivered to flight itinerary option r ∈ R decision variable calculating the arrival time of vehicle k ∈ K at node i ∈ V a large constant number for modeling

θdop , θdcl

Qr Fir ti j ci j xki j ykir aki M

of customer i is shipped by flight itinerary r ∈ R with vehicle k ∈ K. The arrival time of the vehicle k at node j ∈ V is denoted as akj . For the depot, we set akd = θd for any vehicle k ∈ K. Table 1 presents the summary of our notations. The formulation of the ATD-PDP, labeled (MP), is as follows. op



∗ (MP ) zMP = min x,y

  

ci j xki j +

i∈V j∈V \{i}

k∈K





Fir ykir

(1)

i∈C r∈R

Jk (x, y ) =

Subject to



∀i ∈ C,

ykir = 1,

performed by the same vehicle. Constant M is a big number corresponding to arrival times and can be calculated as summation of all the links’ travel times. As mentioned before, if a given flight itinerary option r ∈ R is not useable for a certain customer i ∈ C (e.g., due to not delivering to destination), it is eliminated from further consideration by setting the corresponding yir to zero. For brevity, let Jk (x, y ) denote the objective function for vehicle k and J (x, y ) denote objective for all vehicles as follows.

(2)

 

 ci j xki j + Fir ykir i∈C r∈R  J (x, y ) = Jk (x, y )

i∈V j∈V \{i}

∀k ∈ K

k∈K

k∈K r∈R



xki j ≤ 1,

∀i ∈ V, ∀k ∈ K,

(3)

j∈V \{i}



xki j

i∈V \{ j}



−

xkji

= 0,





 θdcl − trd , Qr , ∀r ∈ R, ∀k ∈ K,

(4)





(5)



aki ≥ ai = θdop + tdi

∀i ∈ C, ∀k ∈ K,

xki j ∀i ∈ C,

∀i ∈ C, ∀r ∈ R, ∀k ∈ K,

∀r ∈ R, ∀k ∈ K,

xkr j ,∀i ∈ C,

∀r ∈ R, ∀k ∈ K,

and, the latest arrival time (a¯ i ) is the latest time that allows a vehicle to pick up the customer load, deliver it to a flight node, and return to the depot before the closing time,

(7)

aki ≤ a¯ i = max min Qr , θdcl − trd − tir

(8) (9)

(10)

j∈V \{r}

ykir , xki j ∈ {0, 1}, aki , akr ≥ 0,

(12)

(6)

j∈V \{i}

ykir ≤

We strengthen the formulation (MP) by network preprocessing and introducing valid inequalities. The first preprocessing step is to tighten the upper and lower bounds on the node arrival times. For a customer node, the earliest arrival time (ai ) is attained via a direct travel from the depot,







r∈R

ykir − 1 M + aki + tir ≤ akr ,

ykir ≤

∀i ∈ V, ∀ j ∈ V \{d, i}, ∀k ∈ K,

∀i ∈ C, ∀k ∈ K,

op , d

akr ≤ min

∀ j ∈ V, ∀k ∈ K,

i∈V \{ j}

xki j − 1 M + aki + ti j ≤ akj ,

aki ≥ θ





3.4. Network Preprocessing and Valid Inequalities

∀i, j ∈ V |i = j, ∀r ∈ R, ∀k ∈ K. (11)

The objective (1) minimizes the total cost of delivery including flight itineraries, service level, and road travel cost. Constraint set (2) ensures that every customer’s load is assigned to a flight itinerary. Constraint set (3) guarantees that each node is visited at most once by each vehicle. Constraint set (4) is the flow conservation at each node for each vehicle. Constraint sets (5) and (6) calculate the arrival time at every node while also preventing subtours. Constraint set (7) prohibits visiting a flight node after the departure time of the flight itinerary while ensuring that the vehicle can also return to the depot before the depot’s closing time. Constraint set (8) guarantees that a customer load pickup precedes its delivery to the selected flight node. Constraint set (9) and (10) ensures that both pickup and delivery of a customer load is



∀i ∈ C, ∀k ∈ K.

(13)

For a flight node, the earliest arrival time (ar ) is attained by the shortest travel from the depot after visiting a customer,

akr ≥ ar = θdop + min (tdi + tir ) i∈C

∀r ∈ R, ∀k ∈ K.

(14)

The latest arrival time to a flight node (a¯ r ) is already included in (MP) as the constraint set (7). The next preprocessing step is determining the lowest value for M in constraints (5) and (8). In particular, we replace M with edge specific Mi j ,

Mi j = a¯ i + ti j

∀i ∈ V, ∀ j ∈ V \{d, i}, ∀k ∈ K.

(15)

The final preprocessing step is the elimination of the inadmissible edges that can never be traversed in a feasible solution. We remove edges (i, d ) ∀i ∈ C since vehicles must return to the depot empty. The edges (d, r ) ∀r ∈ R are eliminated since vehicles leave the depot empty. Lastly, we remove any edge (i, j ) ∀i, j ∈ V if ai + ti j > a¯ j , i.e., the vehicle cannot traverse an edge if it cannot arrive its destination before the latest allowed arrival time. We tighten the constraints set (5) by using the following lifting scheme from Desrochers and Laporte (1991) by taking the reverse arcs into account.









xkji M − ti j − t ji + xki j − 1 M + aki + ti j ≤ akj

∀i, j = i ∈ V \{d}, ∀k ∈ K

(16)

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

In addition, we introduce the following cut set that ensures that a vehicle visits a customer only if it delivers the customer’s load to a flight node.



xkr j =



j∈V \{i}

∀i ∈ C, ∀k ∈ K

ykir

(17)

r∈R

Considering that the vehicles are assumed to be identical, the solution performance of the model can be further improved by introducing symmetry breaking valid cuts. We propose two sets of cuts as follows.

 

 

ci j xki j ≤

i∈V j∈V \{i}

i∈V j∈V \{i}

ci j xki j+1 ∀k ∈ 1 · · · (|K | − 1 )

(18)



2|C |−i τik ≥

i∈C



2|C |−i τik+1

∀k ∈ 1, · · · , (|K | − 1 ),

(19)

i∈C

where τik =

 r∈R

ykir for all i ∈ C and k ∈ K. We set τik to one for the

first customer and first vehicle. The main advantage of these cuts to the previous one is that they do no depend on problem parameters. Our empirical study results demonstrate the notable effectiveness of these cuts in reducing solution times. 4. Solution methodology The (MP) formulation, presented in the previous section, shows that constraint (2) is the only constraints that couple vehicles together. Traditionally, models with similar structures are solved using Lagrangian decomposition. In this section, first, we briefly present the standard Lagrangian Decomposition approach. Next, we discuss the shortcomings of this method and introduce an alternative solution approach, the Successive Subproblem Solution (SSS) method, for solving the ATD-PDP problem using the (MP) formulation. 4.1. Standard Lagrangian decomposition approach The standard Lagrangian Decomposition approach is commonly used for formulations composed of two or more intertwined subproblems that are easier to solve independently through specialized algorithms. This approach has been used for the vehicle routing problems (e.g., Kohl & Madsen, 1997). The (MP) formulation is a candidate for this approach since constraints (2) are the only coupling constraints for vehicles and the rest of the constraints and the objective function is separable by vehicle. Hence, by relaxing constraints (2) through Lagrangian relaxation, (MP) can be decomposed to |K | subproblems, each corresponding to a single vehicle. The Lagrangian relaxation of (MP) with respect to constraints (2) results in the following Lagrangian dual function,



(λ )



=min (x,y )∈

Jk (x, y ) +

k∈K

 i∈C



λi 1 −





ykir

,

(20)

k∈K r∈R

where λ = (λ1 , · · · , λ|C| ) ∈ R|C| is the vector of Lagrangian multipliers associated with constraints (2) and  is the set of feasible solutions. Then, the Lagrangian Dual problem (LD) can be defined as ∗ . follows; indeed, ∗LD is a lower bound on zMP

(LD ) ∗LD = max ((λ ) ). λ

The set of feasible solutions to (MP) can be decomposed into |K | disjoint subsets, i.e.  = 1 × 2 × · · · × |K| , where each k is defined by constraints (3)–(10) for a given k ∈ K. Let us define Lk (λ, x, y ) for each k ∈ K as follows.

Lk (λ, x, y ) = Jk (x, y ) −

(21)

λ i ykir

Consequently, the Lagrangian dual function (λ ) can be presented as follows.

(λ ) =



k (λ ) +



λi ,

i∈C

where k (λ ) = min(x,y)∈k Lk (λ, x, y ) is the Lagrangian function of

the kthth subproblem. To solve (LD), we solve the primal subproblem k (λ ) for each vehicle k at the low-level and update the Lagrangian multipliers at the high-level, e.g., using subgradient optimization (Conejo, Castillo, Minguez, & Garcia-Bertrand, 2006; Fisher, 2004; Geoffrion, 1974). The optimization at both levels is performed iteratively until the dual solution converges. However, since the vehicles are homogeneous, the subproblems k (λ ) are identical for all k ∈ K; i.e. 1 = 2 = · · · = |K| . Hence, all subproblems have the same optimal solution with identical objective values. Accordingly, for solving (LD) at each iteration we need to only solve one subproblem for an arbitrary vehicle k ∈ K and multiply the result by |K | to aquire the lower level solution; that is





∗ LD

= max λ

|K | min Lk (λ, x, y) +



(x,y )∈k



λi

(22)

i∈C

While this solution method reduces the complexity of the problem by decomposing it to subproblems that are easier to solve, the identical subproblems present challenges in the solution process. In particular with discrete decisions, it leads to oscillating dual solutions, affecting the convergence rate. Further, the solutions converged are often primal infeasible and provide a lower bound on ∗ zMP that can be weak. Finally, the primal infeasibility of the solutions requires integration with an exact (heuristic) method such as branch-and-bound method (Lagrangian heuristic) to obtain optimal (good quality) solutions (see for example Kohl & Madsen, 1997). Accordingly, the traditional Lagrangian decomposition approach is not appropriate for solving this problem efficiently. In the next section, we propose an alternative approach to overcome the shortcomings of the Lagrangian decomposition method. 4.2. Successive subproblem solving method We adapt the Successive Subproblem Solving (SSS) method to avoid the challenges associated with the Lagrangian Decomposition method due to the identical subproblems. This approach is introduced by Zhai, Guan, and Cui (2002) to solve the unit commitment problem for scheduling electrical power generators. The SSS approach extends and improves over the Lagrangian Decomposition method by addressing the dual solution oscillation. However, it guarantees neither the primal feasibility nor the quality of achieved solutions. We address these issues in Section 4.2.1 by developing a modified variable target value method for subgradient optimization for the SSS approach. In the SSS, we introduce an absolute penalty term that helps to reduce the oscillation and constraint violations more rapidly. For this purpose, the Lagrangian function is revised to the following augmented form,

1

We are grateful to anonymous reviewer for suggesting the second set of symmetry breaking cuts.

 i∈C r∈R

k∈K

These cuts break symmetry by forcing vehicles to have incremental costs. An alternative set of cuts is as follows1 .

7

Lˆ (ω, λ, x, y ) =

 k∈K

Jk (x, y ) +

 i∈C



λi 1 −





ykir

k∈K r∈R

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR 8

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15





 

+ω ykir ,

1 − i∈C

(23)

k∈K r∈R

where ω > 0 is the penalty parameter. The revised dual problem ˆ ω, λ ) are then expressed as, (PS) and dual function (



∗ PS

ˆ ω, λ ) = max (PS )  (ω ) = max ( λ



min Lˆ (ω, λ, x, y ) .

(x,y )∈

λ

(24)

The ∗PS (ω ) is the optimum dual solution with penalty weight ω. The following theorem establishes that the ∗PS (ω ) is a ∗ . lower bound on the primal optimum solution zMP ˆ ω, λ ) ≤ ∗ (ω ) ≤ z∗ ≤ J (x, y ). Theorem 3. For any ω and λ, ( PS MP ∗ ≤ J (x, y ) and Proof. By definition from (1) and (24), we have zMP ∗ ∗ ∗ ˆ (ω, λ ) ≤  (ω ), respectively. Let (x , y ) be the primal optimum PS

λ∗

solution to problem (MP) and denote the optimum multipliers. The primal optimum solution is feasible, thus, satisfies constraint set (2). Accordingly, we have ∗ zMP

=



∗

∗

Jk (x , y ) +





λ

∗ 1

1−

i∈C

k∈K





yirk∗

k∈K r∈R



 

+ω yirk∗ = Lˆ (ω, λ∗ , x∗ , y∗ ).

1 −

i∈C k∈K r∈R From definition (24), ∗  (ω ) = min Lˆ (ω, λ∗ , x, y) ≤ Lˆ (ω, λ∗ , x∗ , y∗ ) = zMP . ∗ PS

(x,y )∈

 While the introduction of the penalty term empirically alleviates the mentioned oscillation issue, the revised Lagrangian function (23) cannot be decomposed into k subproblems due to the ˆ ω, λ ) with augmentation. Hence, to calculate the subgradient of ( respect to λ, we now need to solve the integrated low-level problem min(x,y)∈ Lˆ (ω, λ, x, y), which is computationally inefficient. The revised Lagrangian function in (23), however, can be reformulated as an additive function. Let us redefine the Lagrangian function for kth vehicle as follows:

Lˆk (ω, λ, x, y ) = Jk (x, y ) −







+ω q i − ykir , ( )

k i∈C

where,

qk ( i ) = 1 −

r∈R

gij = 1 −



ykir

j

∀i ∈ C .

(29)

k∈K r∈R

We first introduce the notation used in the SSS method and then present its algorithmic steps. The SSS algorithm uses ∗PS (ω ) that is often unknown. More discussion on this issue is provided at the end of this section and in Section 4.2.1. Notation: (x j , y j )k : solution of kthth subproblem at iteration j (x j , y j ): solution at iteration j ˆ 0 = {λ ˆ 0 , ∀i ∈ C } λˆ 0 : initial Lagrangian multipliers, i.e., λ i λ j: Lagrangian multipliers at iteration j, i.e., λ j = {λij , ∀i ∈ C} j g: surrogate subgradients at iteration j, i.e., g j = {gij , ∀i ∈ C} δ j: step-size at iteration j j Lagrangian function value at iteration j with penalty ω, Lω : j i.e., Lω = Lˆ (ω, λ j , x j , y j ) β: step-size update parameter, 0 < β < 1 α: initialization factor for Lagrangian multipliers ε: threshold for Lagrangian multiplier convergence criteria, ε>0 εGAP : optimality gap threshold SSS Procedure: Initialization. ˆ 0 , e.g., λ ˆ 0 = 0, solve (LD) using (22) and obI.1. Given λ 0 0 tain (x , y ) I.2. Calculate,

λi ykir

i∈C r∈R

other, thus alleviating the issues associated with identical subproblems. In solving (PS), the SSS method solves the vehicle subproblems one at a time, while calculating the qk (i ) ∀i ∈ C using the solution from other vehicles. The SSS method updates the Lagrangian multipliers after solving any of the subproblems. In the SSS, the Lagrangian multipliers are updated using the surrogate subgradient (SSG) approach introduced by Zhao, Luh, and Wang (1999). The standard subgradient approach requires solving all subproblems to obtain the subgradient direction (Fisher, 2004; Geoffrion, 1974). In the SSG approach, however, the solution to only one of the subproblems is sufficient to obtain a proper surrogate subgradient dij rection. Thus, let gi denote the surrogate subgradient for customer i at any iteration j and be calculated as,

(25)



λ0i = α 1 −



k 0



∀i ∈ C,

yir

(30)

k∈K r∈R

 s∈K s=k

ysir .

(26)

r∈R

It can be verified that the Lagrangian function (23) can be expressed in terms of Lˆk (ω, λ, x, y) and qk (i ) as follows:

Lˆ (ω, λ, x, y ) = Lˆk (ω, λ, x, y ) +



Js (x, y ) +

s∈K s=k



λi qk (i ).

(27)

i∈C

Since (27) is additive, we can now solve problem (PS) in parts, e.g., for each vehicle. The subproblem for vehicle k is then defined as follows:

ˆ k (λ ) = min Lˆk (ω, λ, x, y ) . ( PSk )  (x,y )∈k

(28)

The variable qk (i ) is fixed for the kth subproblem. The qk (i ) links subproblem k to other subproblems by conveying the information about customer i’s assignment to other vehicles. Hence, the solutions of the subproblems are likely to be different from each







0 2

where, 0<α <(∗PS (ω )−Lˆ (ω,0,x0 ,y0 ))/ i∈C 1− k∈K r∈R (ykir )

I.3. Calculate L0ω = Lˆ (ω, λ0 , x0 , y0 ) and update Lagrangian multipliers:

λ1 = λ0 + δ 0 g0 , where 0 < δ 0 = β (∗PS (ω ) − L0ω )/ g0 2 and 0<β <1. Set j=1. Step 1. Subproblem solution: 1.1. For k = 1, 2, . . . , |K |, Repeat: 1.1.a. Solve subproblem (PSk ) in (28) by setting (x j , y j )s = (x j−1 , y j−1 )s for s ∈ K s = k to obtain (x j , y j )k . 1.1.b. If the following improvement condition is satisfied,

Lωj = Lˆ





ω, λ j , x j , y j < Lˆ ω, λ j , x j−1 , y j−1 ,

(31)

where (x j , y j ) = (x j , y j )k ∪ {(x j−1 , y j−1 )s |s ∈ K, s = k}, then go to Step 2; otherwise, continue with the next k. 1.2. Set (x j , y j ) = (x j−1 , y j−1 ). Step 2. Subgradient optimization:

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

2.1. Update Lagrangian multipliers:

λ j+1 = λ j + δ j g j , where 0 < δ j = β (∗PS − Lω )/ g j 2 and 0 < β < 1. Step 3. Check the stopping criteria. ¯ j − Lωj ) ≤ εGAP or λ j+1 − λ j ≤ ε , then go to Step 4; 3.1. If ( PS otherwise set j = j + 1 and return to Step 1. Step 4. Terminate with solution (x j , y j ). j

The SSS method is initialized by solving (LD) to obtain initial solutions to estimate the starting values for the Lagrangian multipliers. The bounds on α in the initialization step ensures that L0ω = Lˆ (ω, λ0 , x0 , y0 ) < ∗PS (ω ). This inequality is important for convergence analysis as explained in the next section. The subproblems in Step 1 are sequentially solved until the improvement condition in (31) is attained. In each subproblem solution, the previous iteration’s solutions are used to calculate qk (i ) ∀i ∈ C. When none of the vehicle k’s subproblem solutions satisfies (31), the previous iteration’s solution is maintained. The multipliers are updated using the surrogate subgradient in Step 2.1. The SSS method terminates when the solution converge to the optimal solution or multipliers converge. Another heuristic stopping criterion such as those based on CPU time or number of iterations can also be adopted, if desired. The SSS method requires ∗PS (ω ). This value, however, is generally unknown in advance and needs to be estimated. A poor underestimation may result in convergence to a primal infeasible solution with large duality gap (see Theorem A.1). In the Lagrangian method, the value used in place of ∗PS (ω ) is an overestimation ∗ , which affects the convergence rate. However, the solutions of zMP converged are either primal infeasible or optimal (Held, Wolfe, & Crowder, 1974). In comparison, the SSS method, using an overestimation of ∗PS (ω ), may converge to a primal feasible but not optimal solution. Hence, the SSS differs from the Lagrangian method, as it may converge to a suboptimal primal feasible solution. The bound estimate of ∗PS (ω ) in the SSS is therefore critical and is affecting both the convergence rate and the solutions converged, i.e., primal feasible or infeasible. We present the bound estimation procedure in Section 4.2.1. We provide the convergence results for SSS method with surrogate subgradient optimization in Appendix A. Further discussion on the convergence, uniqueness of optimal multipliers and other characteristics of the surrogate subgradient algorithm is provided by Chang (2008). 4.2.1. Bound estimation: variable target value method The SSS procedure uses an estimate of ∗PS (ω ) for the surrogate subgradient optimization. Rather than using a static estimate, we dynamically change this estimate in order to obtain a good quality primal feasible solution. Specifically, we modify the variable target value method (VTVM) presented in Lim and Sherali (2006) and incorporate backtracking to improve the target value estimation. Since the SSS method can converge to a primal feasible but suboptimal solution, we integrated a backtracking phase within the VTVM to improve the quality of the feasible solution. We modify the SSS method by replacing ∗PS (ω ) with a j

dynamically adjusted estimate PS (target value). Analogous to j

j

Theorem 3, it can be shown that Lω < PS holds true for each j iteration j. In choosing the estimate PS , the goal is to approxi∗ mate zMP as close as possible. In standard VTVM method, the tarj

get value PS is increased as long as the convergence rate is satisfactory and then decreased towards an optimal solution. In our j adaptation, we increase the target value PS with a controlled rate until we find a primal feasible solution. This primal feasible solution, however, may be a low quality suboptimal solution. Therefore,

9

with a backtracking phase, we revise the latest target value to obtain a better primal feasible solution. Specifically, after encountering with a primal feasible solution, we return to a past iteration where the target value underestimates the current solution’s objective value. Then, the modified SSS repeats the iteration with a smaller step size in an effort to find an improved primal feasible solution. We first provide the notation used in VTVM with backtracking and then present the modified steps of the SSS procedure. We use the recommended value by Lim and Sherali (2006) for η and σ . Next, we briefly discuss the convergence behavior of the SSS with j backtracking. Note that we replace ∗PS (ω ) with PS in the remaining steps of the SSS procedure. Notation for SSS with Backtracking VTVM: j PS : target value at iteration j ¯ j :  upper bound on the optimal solution at iteration j PS

PS : lower bound on the optimal solution value (x∗ , y∗ ) : an optimal solution to (MP) j : accumulated improvements since the last Lagrangian

ε j:

function improvement until the beginning of iteration j acceptance tolerance that the current incumbent value j

σ: ηj : εGAP :

j

Lω is close to the target value PS in iteration j acceptance interval parameter fraction of cumulative improvement that is used to increase the target value in iteration j optimality gap threshold

Modified Steps of the SSS Procedure with Backtracking VTVM: Initialization. Execute Steps I.1, I.2, I.3 of the original SSS procedure, and, j=1 j=1 ¯ PS I.4. Set PS = PS ,  = +∞, η j=1 = 0.35, σ = 0.2,  j=1 =

0, and ε j=1 = σ (PS − Lω ). Step 1. Subproblem Solution & Backtracking: 1.1. For k = 1, 2, . . . , |K |, Repeat: 1.1.a. Solve subproblem (PSk ) in (28) by setting (x j , y j )s = (x j−1 , y j−1 )s for s ∈ K, s = k and obtain ( x j , y j )k . Denote ( x j , y j ) = ( x j , y j )k ∪ {(x j−1 , y j−1 )s |s ∈ K, s = k}. 1.1.b. If (x j , y j ) is primal feasible, then i. Set (x∗ , y∗ ) = (x j , y j ), ¯ j = Lωj = J (x j , y j ), ii. Set  PS iii. Set algorithm parameters, variables, and solutions back to iteration v, i.e., ¯ j = Lωj }, where v = max{l : lPS <  PS iv. Set j := v, v. Set β = β /2 and repeat iteration j with updated Lagrange multipliers λ j = λ j−1 + δ j−1 g j−1 . 1.1.c. If the following improvement condition is satisfied, j=1

Lωj = Lˆ

j=0





ω, λ j , x j , y j < Lˆ ω, λ j , x j−1 , y j−1 ,

where (x j , y j ) = (x j , y j )k ∪ {(x j−1 , y j−1 )s |s ∈ K, s = k}, then go to Step 2; otherwise, continue with the next k. 1.2. Set (x j , y j ) = (x j−1 , y j−1 ). Step 2. Subgradient Optimization & VTVM: j j 2.1. If Lω > PS − ε j , then

j j+1 j ¯ PS 2.1.a. Update the target value PS = min{Lω + ε j + ηj j ,  },

2.1.b. Update the threshold ε j+1 = σ (PS − Lω ), 2.1.c. Reset j = 0, 2.1.d. Update ηj+1 = min{2ηj , 1}, j+1 j ¯ j+1 =  ¯ j , ηj+1 = ηj and ε j+1 = ε j . otherwise set PS = PS ,  PS PS j+1

j

j

j−1

2.2. Update j+1 = j + (Lω − Lω )

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR 10

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

2.3. Update Lagrangian multipliers:

λ

j+1

=λ +δ g , j

j j

where 0 < δ j = β (PS − Lω )/ g j 2 and 0 < β < 1. Step 3. Check the stopping criteria: ¯ j − Lωj ) ≤ εGAP or λ j+1 − λ j ≤ ε , then terminate 3.1. If ( PS with Step 4; otherwise set j = j + 1 and go to Step 1. Step 4. Terminate with (x j , y j ). j

j

j

The SSS with backtracking VTVM initializes the target value PS with an underestimation PS of the dual optimal value, e.g., linear relaxation. When the dual solution is primal feasible, we perform the backtracking phase in Step 1.1b. This backtracking helps to improve the quality of the subsequent feasible solutions by reverting j to an iteration v satisfying vPS < Lω and repeat the iteration j with j

smaller step size. As the Lω converges to the target value such that j j Lω is within ε j threshold of PS , then Step 2.1.a updates the target ¯ j . This value based on the accumulated improvement j and  PS

update guarantees that the dual solution and the target value is separated by at least ε j while ensuring that the target value does not exceed the upper bound. The threshold ε j is updated in Step 2.1.b. Choosing large values for σ increases ε j . With higher ε j values, j

we are more likely to consider that the Lω is close to the target value and thus update the target value more frequently and with larger increments (Step 2.1.a). This can result in poor feasible so¯ j might not have decreased suffilutions as the upper bound  PS ciently. In contrast, lower σ values reduce the convergence rate. The required ranges for acceptance interval parameter and fraction of cumulative improvement are σ ∈ (0, 1/3] and ηj ∈ (0, 1] (Lim & Sherali, 2006). The algorithm terminates and returns the best primal feasible solution when the gap between the best feasible solution and the Lagrangian dual function value falls below the optimality gap threshold (εGAP ). 5. Computational experiments We investigate the computational and solution quality performance of our formulation and the proposed approach for solving the ATD-PDP. All experimental runs are conducted on a PC with 3.40 gigahertz Intel processor and 8 gigabytes RAM. We use CPLEX 12.6.1. In the following section, we first present the result of our study on the effectiveness of our formulation, specifically, the symmetry breaking cuts (18) and (19). Next, we report on the computational results of both variants of the SSS method, namely SSS with backtracking VTVM (SSS-B-VTVM) and VTVM base SSS without backtracking (SSS-VTVM). Since the ATD-PDP is a new problem, no benchmark instances are available. We generated a set of test problems varying from small to large problem scenarios (see Appendix B). In generating the data sets, we adhered to the development procedure described in Solomon (1987). The problem scenarios have one depot and one or two airports each with three flight itinerary options for each individual customer; the third option represents the recourse flight itinerary option. For a problem scenario with n = |C | customers and m = |H | airports, we first generate (1 + m + n) locations from a uniform distribution over a square area bounded by [0, 10(1 + m + n)] × [0, 10(1 + m + n)]. Next, we randomly label the nodes as depot, airports and customers to avoid any association between the location and the identity of a node. The travel times between nodes are calculated as the Euclidean distances. The travel cost between two nodes is set equal to their travel time. For each airport h, the departure times Qrh of flights are independent and identically distributed according to a uniform distribution U[ϕ /|K |, θ ] where |K | is the number of available vehicles; ϕ

is the heuristic solution to a TSP problem consisting of the depot (origin), all customers and the selected airport (destination). The heuristic solution is obtained through the greedy next best routing heuristic. The cost of flight itinerary options, Firh , are independent and identically distributed according to a uniform distribution U[10 0, 60 0]. The flight itinerary options are sorted from cheapest to most expensive and assigned to the flight itineraries based on the starting times such that cheaper itineraries start earlier. 5.1. Model formulation evaluation For our first study, we generated experiments with only 7 customers and solved them to optimality. As it will be shown later in this section, CPLEX is often not capable of solving larger problems to optimality. For each experiment, we generated 20 independent instances and solved them “without” and “with” the symmetry breaking cuts (18) and (19). We evaluated the effectiveness of the cuts based on the solution time of the instances. The solution time was restricted to 10 0 0 seconds. For all the instances, the objective values of the solutions “with” and “without” the cuts were the same; thus, we limited the evaluation of the cuts performance to the solution times. Table 2 reports the average, maximum and minimum ratio of the solution times for instances with 7 customers, one or two airports, and 3 to 5 vehicles. It also reports the average solution times. The results in Table 2 demonstrate that the symmetry breaking cuts reduce the solution times significantly. As problem instances get more complex, the effectiveness of the cuts tends to increase. This increase in the effectiveness of cuts is particularly more notable as the number of vehicles increases. Comparing cuts (18) with cuts (19), cuts (19) seem to be more effective in improving the computational performance. However, among the cases that the time saving was significant, only for the instance with 2 airports and 5 vehicles the difference between the cuts were statistically significant. As will be presented in the next section, CPLEX often fails to produce the optimal solution for more complex problems. Our empirical study shows that while cuts (19) seem to reduce solution times in general, in those cases where CPLEX is unable to find the optimum solutions, cuts (18) result in better lower bound than those of (19). Accordingly, we use cuts (19) in our computational experiments to solve the problems with CPLEX due to its overall superiority; however, for cases that CPLEX fails to obtain the optimum solution we use the best of the lower bounds produce by cuts (18) or (19). Indeed, these cuts are irrelevant for the decomposed problem; i.e., subproblems with a single vehicle. 5.2. Solution method evaluation In the second part of the computational study, we conducted experiments using 7, 10 and 15 customer cases. For each scenario, we generated 20 independent instances and solved them using CPLEX, SSS-VTVM, and SSS-B-VTVM. Since there is no prior work on ATD-PDP, we compare the proposed methods’ performances against the CPLEX solutions and lower bounds of the model (MP). We restricted the solution times to 300 seconds for all methods and reported the best feasible solution attained within the time limit for each instance. In total, we solved 360 problem instances. In running CPLEX, we used the default settings except that we turned on the feasibility pump of the mixed integer programming with an emphasis on solution feasibility. The ideal scenario is to compare the methods’ best feasible solution with the optimum solution. However, as can be seen from the results below, optimum solutions were only attainable for small problems. Thus, to evaluate the quality of a solution we

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

11

Table 2 The performance of symmetry breaking cuts.

|C |

|H |

7

1

2

|K |

3 4 5 3 4 5

Ratio of Solution Times of “Without” to

Solution Time Averages (seconds)

“With” Cuts (18)

“With” Cuts (19)

“Without” Cuts

“With” Cuts (18)

“With” Cuts (19)

Avg

Max

Min

Avg

Max

Min

0.9 0.8 0.9 2.3 2.9 3.3

2.6 2.3 2.8 4.7 7.9 8.1

0.3 0.2 0.2 0.4 0.3 0.1

1.3 1.4 1.7 2.3 5.7 7.8

5.1 5.5 9.2 4.8 13.6 19.2

0.5 0.2 0.4 0.4 0.5 0.6

2.5 4.1 5.5 95.3 215.2 123.9

2.2 3.7 3.6 37.8 77.3 36.0

1.4 2.4 1.4 37.3 40.2 13.7

Table 3 Comparative performance of LBl against optimum solution and LBc .

|C | 7

|H | 1

2

10

1

2

15

1

2

|K | 3 4 5 3 4 5 3 4 5 3 4 5 3 4 5 3 4 5

Opt

100 100 100 100 100 100 95 90 80 35.0 25.0 30.0 20 25 15 5 0 0

Gap between LBl and optimum solution (%) Avg.

max

min

1.1 3.9 3.6 2.1 3.7 2.6 1.2 1.9 1.9 4.7 4.7 2.4 0.1 0.5 1.8 2.3 -∗∗ -∗∗

9.2 15.0 16.3 11.7 22.0 20.2 8.3 18.7 22.8 11.9 11.6 6.9 0.3 2.9 7.2 2.3 -∗∗ -∗∗

0.1 0.6 0.0∗ 0.1 0.1 0.0∗ 0.0∗ 0.0∗ 0.0∗ 0.0∗ 1.3 0.0∗ 0.0∗ 0.0∗ 0.0∗ 2.3 -∗∗ -∗∗

LBl > LBc (%), if optimum is unknown

–∗∗∗ –∗∗∗ –∗∗∗ –∗∗∗ –∗∗∗ –∗∗∗ 100 0 0 15.4 73.3 71.4 12.5 33.3 41.2 55.0 75.0 70.0

Average LBl time (sec)

86 125 160 293 295 279 476 462 439 569 597 559 593 596 565 614 625 623

∗

These values are small positive numbers greater than zero. Indicates that the corresponding performance for the problem case is not available since no optimal solution is available from CPLEX for any of the problem instances. ∗∗∗ Indicates that the corresponding performance for the problem case is not available since the optimal solution is available from CPLEX for all the problem instances. ∗∗

used the best lower bound as the alternative benchmark. Obtaining a high quality lower bound for this problem is proven to be a challenge. The lower bound achieved from linear programming relaxation tends to have a large gap and is not useful. Accordingly, we obtained two lower bounds, i.e., (1) the best lower bound produced by CPLEX on the integrated problem and (2) the best lower bound attained by using standard Lagrangian relaxation on the decomposed problem as described in Section 4.1. We denote these lower bounds as LBc and LBl , respectively. Indeed, if CPLEX succeeds to solve a problem instance to optimality, the LBc is the same as the optimum solution; otherwise, there will be some gap between LBc and the optimum solution. As discussed in Section 4.1, Lagrangian decomposition does not provide a feasible solution for this problem in most practical cases, thus, LBl is expected to have a gap. We also note that in our experimental problem instances, Lagrangian decomposition never produced a feasible solution. Although LBl is obtained by solving the relaxed subproblems, due to the complexity of the subproblems, it is still a challenge to obtain a high quality lower bound. This is because the subproblems are essentially single vehicle versions of the original problem with the added component of prize collection for visited customers where visiting all customers is not required. We record the best lower bound attained within 600 seconds as LBl ; that is, no new iteration is initiated after time limit is exceeded but ongoing iterations are not terminated. We increased the solution time window for the lower bound calculation since

the Lagrangian decomposition tend to have a slow convergence rate, as discussed in Section 4.1, and requires a longer time to produce a quality lower bound. Table 3 presents the performance of LBl against the optimum solution when available and against the LBc when the optimum solution is not available. In this table, we report the percentage of instances for which the optimum solution is attained (labeled opt). For these instances, we present the average, maximum and minimum gap between the LBl and the known optimum solution. For the remaining instances for which the optimum solutions could not be attained, we report the percentage of times that LBl performed better than LBc . As presented in Table 3, the average gap of LBl is less than 5%. However, the gap can be as much as 23% in some instances. On the other hand, for instances for which the optimum solution cannot be attained, it is impossible to estimate the gap. However, we have observed that LBl tends to be a better choice than LBc . This is particularly applicable to more complex problem instances. In the remainder of this section, we refer to the best of LBc and LBl , simply as the lower bound. In Table 4, we present the results of our computational experiments for three methods: CPLEX, SSS-VTVM, and SSS-B-VTVM. In our experiments, SSS-VTVM and SSS-B-VTVM always reached a primal feasible solution within the time limit. However, CPLEX, particularly for more complex problems, often failed to generate a feasible solution. Thus, for CPLEX, we report the percentage of the cases that CPLEX reached a primal feasible solution (labeled Fsbl)

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR 12

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15 Table 4 The comparative performance of CPLEX, SSS-VTVM, and SSS-B-VTVM.

|C | 7

|H | 1

2

10

1

2

15

1

2

|K | 3 4 5 3 4 5 3 4 5 3 4 5 3 4 5 3 4 5

CPLEX

SSS-B-VTVM

SSS-VTVM

Avg. solution time (seconds)

Fsbl

Opt

AOG

Opt

Avg

Max

Opt

Avg

Max

CPLEX

SSS-B-VTVM

SSS-VTVM

100 100 100 100 100 100 100 100 100 90 90 90 100 95 95 25 35 35

100 100 100 100 100 100 95 90 80 35 25 30 20 25 15 0 0 0

0.2 0.1 0.2 0.2 0.2 0.0 0.2 0.2 0.2 0.3 0.0 0.0 0.1 0.7 0.1 -∗ -∗ -∗

65.0 70.0 70.0 75.0 80.0 95.0 68.4 61.1 75.0 85.7 100 100 75.0 40.0 66.7 -∗ -∗ -∗

0.2 0.1 0.2 0.2 0.2 0.0 0.3 0.2 0.3 2.8 4.7 4.2 4.1 5.1 5.6 9.4 11.0 11.3

1.1 0.6 1.5 1.5 1.5 0.1 3.3 1.8 3.2 11.6 11.5 11.6 8.2 15.9 13.4 15.4 22.7 17.9

35.0 40.0 30.0 50.0 60.0 85.0 26.3 27.8 62.5 71.4 80.0 66.7 0.0 40.0 0.0 -∗ -∗ -∗

2.5 1.6 2.4 0.7 1.1 0.4 2.3 3.7 2.2 4.0 5.8 6.2 7.6 7.7 8.8 10.8 12.5 13.2

10.9 7.1 13.3 5.3 8.1 7.3 11.8 13.7 10.8 11.6 13.9 21.7 14.6 17.8 18.0 17.8 22.8 29.6

1 2 1 18 22 11 68 64 102 231 271 252 263 242 258 300 300 300

5 7 6 39 27 16 22 27 12 70 103 81 63 62 54 136 139 95

1 2 1 13 12 6 7 6 4 35 51 30 12 16 19 63 62 29

∗ Indicates that the corresponding performance for the problem case is not available since no optimal solution is available from CPLEX for any problem instances.

and the percentage of cases for which an optimal solution is obtained (labeled Opt). For SSS-VTVM and SSS-B-VTVM, we present the average (labeled Avg) and maximum (labeled Max) gap against the lower bound. Furthermore, we report the percentage of cases that SSS-VTVM or SSS-B-VTVM achieved an optimal solution. Note that we can conclude optimality for SSS-VTVM and SSS-B-VTVM solutions only if an optimal solution is available through CPLEX. Consequently, in calculating this measure, we ignore the cases where CPLEX failed to obtain an optimal solution. To provide a better perspective on the optimality gap of SSS-B-VTVM solutions, we reported the average optimality gap for instances for which an optimum solution is identified (labeled AOG). Table 4 also presents the average solution times for each of the three methods. The results show that CPLEX is able to solve all the instances of the problem cases with 7 customers and one airport to optimality. However, as the number of customers increases, particularly in combination with an increasing number of airports, the problems get more complex for CPLEX to solve. CPLEX’S ability to achieve optimal solution fades quickly, and for problems with 10 and 15 customers, CPLEX often fails to produce even a feasible solution. The results clearly show that using CPLEX is not a practical approach for solving this problem even for mid-size problems. As evident in Table 4, both SSS-VTVM and SSS-B-VTVM are able to produce high-quality solutions within a notably short time. Among all the cases where an optimum solution is identified, SSSB-VTVM managed to find the same solution in more than 65, 61 and 40 percent of instances for scenarios with 7, 10 and 15 customers, respectively. Moreover, the average optimality gap (AOG) for SSS-B-VTVM for all the cases is less than 0.7 percent. In other words, for any problem instance that we actually have the optimum solution to compare with, SSS-B-VTVM produces a feasible solution that is either the optimum or very close to it. As mentioned earlier, CPLEX failed to produce optimum solutions for many instances. Thus, in absence of optimum solutions, we have to rely on an indirect measure of solution quality by comparing SSS-B-VTVM solutions against best lower bounds. As presented in Table 4, the average gap against the lower bound for all the instances with 7 customers and the instances with 10 customers and one airport is less than 0.3 percent. The worst-case scenario for this group of instances has a gap of less than 3.3 percent. For more complex problems, the average gap increases to about 11 percent with a worst-case scenario of 22.7 percent gap

against the lower bound. However, for understanding these numbers we need to consider that the lower bounds are just estimations of the optimum solutions and as mentioned earlier maybe far away from it. To investigate this possibility further, we conducted an experiment. The worst gap for SSS-B-VTVM solutions against lower bound belongs to a case with 15 customers, 2 airports, and 4 vehicles. The gap was 22.7 percent. We reevaluated this case by allowing extra time for the lower bound calculation. Accordingly, we allowed CPLEX and the standard Lagrangian algorithm to run for 10,0 0 0 seconds to produce potentially better lower bounds. Both methods were terminated due to reaching the time limit without convergence to an optimal solution. At the end of the experiment, the CPLEX’s lower bound did not improve notably. However, LBl improved significantly. Consequently, SSS-B-VTVM solution gap against the lower bound for this instance dropped from 22.7 percent to 4.1 percent. This result confirms that the gap against the lower bound can only be used as an overestimation of the optimality gap for a solution method and the actual optimality gap may be much smaller than these estimations. A similar analogy can be drawn for SSS-VTVM. Overall, the quality of SSS-VTVM solutions is less than those obtained by SSSB-VTVM, as expected. However, the solution times are notably shorter. The average solution time for SSS-VTVM is about one minute or less, compared to 2.4 minutes or less for SSS-B-VTVM. Nevertheless, the improvement in solution quality due to the backtracking phase of SSS-B-VTVM seems to be worth the extra solution time. Overall, the results of the experimental study confirm the performance of our proposed solution methods in terms of both solution time and solution quality. In assessing the solution times, it should be noted that over 90% of the solution times are spent in solving the subproblems. Upon the availability of an efficient method to solve the subproblems, it can be expected that the performance of SSS-VTVM and SSS-B-VTVM methods can be further improved. Possibly, the most noteworthy aspect of the proposed methods is their ability to generate an implementable good quality feasible solution within a short time that is helpful as a decision support tool for practical purposes 6. Conclusion We study a freight forwarder’s operational implementation of AAAP in a MAR for air cargo transportation. The forwarder’s AAAP

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

implementation involves the task of selecting flight itinerary options for a given set of heterogeneous air cargo customers, picking up their loads via a fleet of vehicles and then delivering them to the chosen airports in the region. The goal is to minimize the total cost of air and road transportation by simultaneously selecting the air cargo flight itinerary and scheduling pickup and delivery of multiple customer loads. We formulated a novel model (ATPPDP) which extends the existing PDP models to address the case where the delivery cost is both destination and time dependent. This model is further strengthened by preprocessing steps and special cuts. To overcome the computational complexities, we adapted the SSS to form an efficient solution algorithm. The proposed algorithm overcomes the challenges associated with identical subproblems and iteratively solves the ATP-PDP in parts. Moreover, we developed a modified variable target value methodology for the SSS to enhance the quality of the solutions. The integrated method, SSS-B-VTVM, converges to a primal feasible solution and the solution quality can be controlled by trading-off the quality with computational performance. We conducted an experimental study to assess the optimality gap and CPU time performance of SSS-B-VTVM and compared it with those of CPLEX and lower bounds. The results show that SSSB-VTVM yields near-optimal primal feasible solutions. Further, SSSB-VTVM is able to achieve this performance on average in less than one minute. This is in contrast to CPLEX that failed in most cases to obtain a feasible solution within 300 seconds. In this paper, we introduced a new class of PDPs as ATP-PDP. Future research can focus on expanding this class and obtaining an effective heuristic algorithm for solving the ATP-PDP. One avenue of future research is to exploit the special structure of the subproblems for individual vehicles to develop a more efficient solution method for the single vehicle problems. Another extension of the proposed modeling and algorithmic approach is to consider capacitated heterogeneous fleet where the variations amongst trucks would affect the assignment and routing decision. Acknowledgments

Appendix A Theorem A.1. (Solution Bounding) For a given ω, at each iteration j i, Lω < ∗PS (ω ). Proof. For j = 0, the condition on α suffices. In the case of j ≥ 1, from (31) we have,

Lωj = Lˆ Lˆ



The following lemma states that the search direction of the Lagrangian multipliers in any iteration is always a proper direction, i.e., (λ∗ − λ j )g j > 0. Lemma A.1. (Direction). Let λ∗ be the optimal multiplier vector, then ∗PS (ω ) − Lωj ≤ (λ∗ − λ j )g j , ∀ j. Proof. Based on (23) and (24), we have



ˆ ω, λ∗ ) = Lˆ (ω, λ∗ , x∗ , y∗ ) ≤ Lˆ ω, λ∗ , x j , y j ∗PS (ω ) = (

= Lωj + Lˆ ω, λ∗ , x j , y j − Lωj

= Lωj + λ∗ − λ j g j . The last step follows from the definition of g j in (29) and Lagrangian function Lˆ (ω, λ, x, y) in (23). From Theorem A.1, we have ∗PS (ω ) − Lωj > 0, thus the theorem’s result follows.  The convergence of the Lagrangian multipliers is established by the following theorem. Theorem A.2. (Convergence) In the SSS algorithm, the Lagrangian multipliers are converging; i.e.,

 ∗    λ − λ j+1 2 < λ∗ − λ j 2 ∀ j, where λ∗ is the optimal multiplier vector. Proof. From (32) we have

 ∗    λ − λ j+1 2 = λ∗ − λ j − δ j g j 2  2 2  2

= λ∗ − λ j  + δ j g j  − 2δ j λ∗ − λ j g j . Using result from Lemma A.1, we have,

 ∗     

λ − λ j+1 2 ≤ λ∗ − λ j 2 + δ j 2 g j 2 − 2δ j ∗PS (ω ) − Lωj . Then, from the definition of δ j in Step 2 of SSS procedure,

We are thankful to the anonymous reviewers for their invaluable comments and suggestions, which significantly improved the clarity, technical contents, and presentation of this paper. This work was supported by funding grant DTRT06-G-0039 from the US Department of Transportation through the University Transportation Center at University of Toledo (UT-UTC).



13

j



ω, λ j , x j , y ≤ Lˆ ω, λ j , x j−1 , y

j−1

= Lω

+δ

j−1

 j−1 2 g  .

Increasing the penalty parameter improves the quality of the solution converged as established by the following theorem. Theorem A.3. For any two penalty weight ω1 and ω2 , where 0 < ω1 < ω2 ,

ˆ ω2 , λ ) ≤ ∗ (ω2 ) ≤ z∗ . ˆ ω1 , λ ) ≤ ( ( PS MP Proof. From (23), we have L j ω − L j ω = (ω2 − ω1 )| j

Further,

i∈C

and using the result of Theorem A.1, we obtain λ∗ − λ j+1 2 ≤

λ∗ − λ j 2 . 

j

2



1

i∈C

gi | ≥ 0. j

Thus, Lω2 ≥ Lω1 . Subsequently from (24), we have,

.



ω, λ j , x j−1 , y j−1 = Lˆ ω, λ j−1 , x j−1 , y j−1



+Lˆ ω, λ j , x j−1 , y j−1 − Lˆ ω, λ j−1 , x j−1 , y j−1    j

j−1 j−1 k = Lω + 1− yir λi − λi j−1

 ∗   

λ − λ j+1 2 ≤ λ∗ − λ j 2 − δ j ∗PS (ω ) − Lωj ,

minLωj 2 ≥ minLωj 1 , ˆ ω2 , λ ) ≥ max( ˆ ω1 , λ ), max(

∗PS (ω2 ) ≥ ∗PS (ω1 ). ∗ . From Theorem A.2, we already have ∗PS (ω2 ) ≤ zMP



k∈K r∈R

From the definition of δ j in Step 3 of SSS procedure we have, j j−1 j−1 j j−1 Lω ≤ Lω + β (∗PS (ω ) − Lω ). Since β < 1, we obtain, Lω < Lω + ∗PS (ω ) − Lωj−1 ≤ ∗PS (ω ). 

Appendix B Problem instances used in Section 5 along with their descriptions are available online as supplementary material with this article.

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

ARTICLE IN PRESS

JID: EOR 14

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15

Appendix C

flights are useable for a given load’s requirements. As explained in Section 3, we can easily adopt these restrictions in our model.

In this appendix, we provide a realistic example to illustrate the air cargo transportation under the AAAP. For this example, we use actual data for a scenario with cargo shipping out of southern California MAR. We extracted the data from the public database of the Bureau of Transportation Statistics (BTS) accessible at http://www.transtats.bts.gov. We present this example with only one cargo carrier, Southwest Airlines (WN). Southwest airlines flies out of four airports in southern California MAR. They are: Bob Hope Airport (BUR), Los Angeles International Airport (LAX), Ontario International Airport (ONT), and John Wayne Airport (SNA). These airports are in close vicinity of each other (50 miles or less) and constitute a MAR. It is common for air cargo to be shipped through connecting flights to its destination. For this example, however, we only consider direct flights. From the south California MAR, 38 destinations (single airport or MARs) are accessible with direct flights by Southwest Airlines. Among these destinations, 18 destinations (47%) are accessible from more than one airport in the origin MAR. Furthermore, five destinations are accessible from all the airports in the MAR. This latter group of destinations accounts for 60% of freight traffic (by weight) out of the total freight transported out of the MAR by Southwest airlines to any domestic destination (for November 2014). There are multiple flights from each origin airports in the MAR to each destination (as many as 34 per day). Table C.1 presents a set of flights departing from the origin MAR to one of the common destinations for a given weekday (departure time between 13:00 and 18:00; November 19th, 2014). In this table, cut-off times are the latest cargo acceptance time ahead of the scheduled departure time. For instance, the latest time to deliver a load to be carried to LAS by the 15:10 flight out of LAX is 14:10. As presented in Table C.1, there are 10 flight itinerary options (in the chosen time window) to ship cargo to LAS. Each flight itinerary option has a drop-off cut-off time and will result in a different arrival time at the destination. For a given customer, the cost of each flight itinerary option is different due to the arrival time and itinerary shipping rates. The cost of a single flight itinerary option is also generally different for two customers due to factors such as their individual chargeable weights, cargo class, transportation requirements, and cost of arrival delays. In this example, we only consider direct flights. If we include destinations with connections, we could have a situation that two customers are using the same flight leg out of an origin airport, thus subject to the same cut-off times, but continue to their individual destinations through different consequent flights, thus subject to different flight itinerary costs. Similarly, not all the destinations are accessible from each origin airport in a MAR or not all

Table C.1 An example of daily flight schedule from south California MAR to a common destination. Origin

Destination

Scheduled departure from origin

Scheduled arrival at destination

BUR LAX SNA ONT LAX BUR LAX SNA BUR ONT

LAS LAS LAS LAS LAS LAS LAS LAS LAS LAS

13:00 13:40 13:50 15:00 15:10 15:30 16:35 16:40 17:00 18:00

14:00 14:45 14:50 16:00 16:10 16:30 17:40 17:45 18:05 19:00

Cut-off time 30 60 45 30 60 30 60 45 30 30

Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ejor.2017.04.033. References Airbus (2015). Global Market Forecast 2015-2034 , Airbus (Ed.), Blagnac Cedex, France. Anily, S., & Bramel, J. (1999). Approximation algorithms for the capacitated traveling salesman problem with pickups and deliveries. Naval Research Logistics, 46(6), 654–670. Anily, S., & Hassin, R. (1992). The swapping problem. Networks, 22(4), 419–433. Azadian, F., Murat, A. E., & Chinnam, R. B. (2012). Dynamic routing of time-sensitive air cargo using real-time information. Transportation Research Part E: Logistics and Transportation Review, 48(1), 355–372. Berbeglia, G., Cordeau, J.-F., Gribkovskaia, I., & Laporte, G. (2007). Static pickup and delivery problems: A classification scheme and survey. TOP, 15(1), 1–31. Berbeglia, G., Cordeau, J.-F., & Laporte, G. (2010). Dynamic pickup and delivery problems. European Journal of Operational Research, 202(1), 8–15. Boeing (2015). Current Market Outlook 2015-2034 , Boeing (Ed.), Seattle, WA. Chalasani, P., & Motwani, R. (1999). Approximating capacitated routing and delivery problems. SIAM Journal on Computing, 28(6), 2133. Chang, T. S. (2008). Comments on “Surrogate gradient algorithm for Lagrangian relaxation”. Journal of Optimization Theory and Applications, 137(3), 691–697. Conejo, A., Castillo, E., Minguez, R., & Garcia-Bertrand, R. (2006). Decomposition techniques in mathematical programming engineering and science applications. Berlin; Heidelberg; New York: Springer. Cordeau, J.-F., & Laporte, G. (2007). The dial-a-ride problem: Models and algorithms. Annals of Operations Research, 153(1), 29–46. Desrochers, M., & Laporte, G. (1991). Improvements and extensions to the Miller– Tucker-Zemlin subtour elimination constraints. Operations Research Letters, 10(1), 27–36. Diana, M., & Dessouky, M. M. (2004). A new regret insertion heuristic for solving large-scale dial-a-ride problems with time windows. Transportation Research Part B: Methodological, 38(6), 539–557. Figliozzi, M. A. (2007). Analysis of the efficiency of urban commercial vehicle tours: Data collection, methodology, and policy implications. Transportation Research Part B: Methodological, 41(9), 1014–1032. Fisher, M. L. (2004). The Lagrangian relaxation method for solving integer programming problems. Management Science, 50(12), 1861–1871. Frederickson, G. (1978). Approximation algorithms for some routing problems. SIAM Journal on Computing, 7(2), 178. Geoffrion, A. M. (1974). Lagrangian relaxation for integer programming. In M. L. Balinski (Ed.), Approaches to Integer Programming (pp. 82–114). Berlin Heidelberg: Springer. Hall, R. W. (2002). Alternative Access and Locations for Air Cargo. METRANS Transportation Center, University of Southern California. Held, M., Wolfe, P., & Crowder, H. P. (1974). Validation of subgradient optimization. Mathematical Programming, 6(1), 62–88. Hellermann, R. (2006). Capacity options for revenue management theory and applications in the air cargo industry. Berlin; Heidelberg; New York: Springer. Hernández-Pérez, H., & Salazar-González, J.-J. (2004a). A branch-and-cut algorithm for a traveling salesman problem with pickup and delivery. Discrete Applied Mathematics, 145(1), 126–139. Hernández-Pérez, H., & Salazar-González, J.-J. (2004b). Heuristics for the one-commodity pickup-and-delivery traveling salesman problem. Transportation Science, 38(2), 245–255. Hernández-Pérez, H., & Salazar-González, J.-J. (2007). The one-commodity pickup-and-delivery traveling salesman problem: Inequalities and algorithms. Networks, 50(4), 258–272. Hernández-Pérez, H., & Salazar-González, J.-J. (2009). The multi-commodity one– to-one pickup-and-delivery traveling salesman problem. European Journal of Operational Research, 196(3), 987–995. Hernández-Pérez, H., & Salazar-González, J. J. (2014). The multi-commodity pickup-and-delivery traveling salesman problem. Networks, 63(1), 46–59. Kirchler, D., & Calvo, R. W. (2013). A granular Tabu search algorithm for the dial-a-ride problem. Transportation Research Part B: Methodological, 56, 120–135. Kohl, N., & Madsen, O. B. G. (1997). An optimization algorithm for the vehicle routing problem with time windows based on Lagrangian relaxation. Operations Research, 45(3), 395–406. Laporte, G. (1992). The vehicle routing problem: An overview of exact and approximate algorithms. European Journal of Operational Research, 59(3), 345–358. Laporte, G. (2009). Fifty years of vehicle routing. Transportation Science, 43(4), 408–416. Lim, C., & Sherali, H. (2006). Convergence and computational analyses for some variable target value and subgradient deflection methods. Computational Optimization and Applications, 34(3), 409–428. Loo, B. P. Y. (2008). Passengers’ airport choice within multi-airport regions (MARs): some insights from a stated preference survey at the Hong Kong international airport. Journal of Transport Geography, 16(2), 117–125.

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033

JID: EOR

ARTICLE IN PRESS

[m5G;May 2, 2017;13:24]

F. Azadian et al. / European Journal of Operational Research 000 (2017) 1–15 Paquette, J., Cordeau, J. F., Laporte, G., & Pascoal, M. M. (2013). Combining multicriteria analysis and Tabu search for dial-a-ride problems. Transportation Research Part B: Methodological, 52, 1–16. Parragh, S., Doerner, K., & Hartl, R. (2008a). A survey on pickup and delivery problems (Part I: Transportation between customers and depot). Journal für Betriebswirtschaft, 58(1), 21–51. Parragh, S., Doerner, K., & Hartl, R. (2008b). A survey on pickup and delivery problems (Part II: transportation between pickup and delivery locations). Journal für Betriebswirtschaft, 58(2), 81–117. Savelsbergh, M. W. P. (1985). Local search in routing problems with time windows. Annals of Operations Research, 4, 285–305.

15

Solomon, M. M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations Research, 35(2), 254–265. Toth, P., & Vigo, D. (2001). The Vehicle Routing Problem. Philadelphia: Society for Industrial and Applied Mathematics. Zhai, Q., Guan, X., & Cui, J. (2002). Unit commitment with identical units successive subproblem solving method based on Lagrangian relaxation. IEEE Transactions on Power Systems, 17(4), 1250–1257. Zhao, X., Luh, P. B., & Wang, J. (1999). Surrogate Gradient Algorithm for Lagrangian Relaxation. Journal of Optimization Theory and Applications, 100(3), 699–712.

Please cite this article as: F. Azadian et al., An unpaired pickup and delivery problem with time dependent assignment costs: Application in air cargo transportation, European Journal of Operational Research (2017), http://dx.doi.org/10.1016/j.ejor.2017.04.033