Discrete intermodal freight transportation network design with route choice behavior of intermodal operators

Transportation Research Part B 95 (2017) 76–104

Contents lists available at ScienceDirect

Transportation Research Part B journal homepage: www.elsevier.com/locate/trb

Discrete intermodal freight transportation network design with route choice behavior of intermodal operators Xinchang Wang a, Qiang Meng b,∗ a b

Department of Marketing, Quantitative Analysis, and Business Law, Mississippi State University, MS 39762, USA Department of Civil and Environmental Engineering, National University of Singapore, Singapore 117576, Singapore

a r t i c l e

i n f o

Article history: Received 20 July 2016 Revised 25 September 2016 Accepted 2 November 2016

Keywords: Freight transportation Intermodal network Economies of scale Congestion effects Route choice model Global optimization algorithm

a b s t r a c t We consider a discrete intermodal network design problem for freight transportation, in which the network planner needs to determine whether or not to build up or expand a link to minimize the total operating cost of carriers and hub operators under a general route choice model of intermodal operators. We formulate the problem as a mixed-integer nonlinear and non-convex program that involves congestion effects, piecewise linear cost functions, and a fixed-point constraint. We develop a series of relaxed and equivalent models to reduce the hardness of the problem and provide theoretical results to show the equivalences. We present two solution methods to solve the problem with one returning heuristic solutions and the other generating a globally optimal solution. We offer two numerical experiments to test the two solution algorithms and also shed light on their performance comparisons. © 2016 Published by Elsevier Ltd.

1. Introduction In this paper, we look into the discrete intermodal transportation network design (ITND) problem, in which there exists a network planner who attempts to re-design an existing intermodal freight transportation network, while taking care of the benefits of other stakeholders involved, such as carriers and hub operators, and taking into account the route choice behavior of intermodal operators. It is worth noting that, we do not restrict the applications of our study to the setting where a network has to pre-exist. As will be seen later in Section 4, our study can be applied to a special case where there is no existing network and the network planner needs to build up a network from scratch. We will use “cargo” and “freight” interchangeably in this paper. The network planner could represent a regional, national, or international government division that has made a plan to improve/design an intermodal freight transportation network to foster intermodal cargo flows within its governed area. One example is a port city that intends to develop a dedicated hinterland network to help cargo flows access the port in a smoother and cheaper fashion (Wang et al., 2016). Another example is the Economic and Social Commission for Asia and the Pacific (ESCAP) of the United Nations which has advocated building an integrated transport network, including road, rail, dry ports, and their efficient linkages to seaports and airports, within its member countries (ESCAP, 2007). Carriers refer to land transport companies or ocean shipping companies that provide cargo transport services using vehicles, crew, and equipment. ∗

Corresponding author. E-mail addresses: [email protected] (X. Wang), [email protected] (Q. Meng).

http://dx.doi.org/10.1016/j.trb.2016.11.001 0191-2615/© 2016 Published by Elsevier Ltd.

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

77

Hub operators could be thought of as operating companies that provide cargo transfer services at intermodal hubs, such as dry ports, seaports, and border crossing terminals between two adjacent countries, at which cargo needs to be transshipped from one country’s vehicles to the vehicles of the other one. We herein use intermodal operators to represent all types of entities who make mode and route selections for cargo shipments that need to be transported through the intermodal network and coordinate the entire intermodal transportation processes for cargo flows between origin-destination (O–D) pairs (Meng and Wang, 2011). Examples of intermodal operators include freight forwarders, the third-party logistics companies, and shippers having their own transport departments. Different stakeholders intervene in the ITND problem with quite different characteristics and behaviors. We now describe how each stakeholder participates in the play. Carriers and hub operators. Carriers and hub operators are end users of network infrastructures such as vehicles, road lanes, and rail tracks. They could bear a fixed installation cost even when they handle no cargo flows. Moreover, the operations of carriers and hub operators are often believed to exhibit the so-called economies of scale, which refer to the phenomenon that the unit cargo handling cost declines with the growth of the total volume of cargo transported/transferred (Kimms, 2006). We will, therefore, formulate the cost functions for carriers and hub operators to be piecewise linear functions with decreasing segment slopes to reflect the economies of scale and with a constant term to model the fixed installation cost. The network planner. One of the prominent features that distinguish the ITND problem from the traditional unimodal transportation network design problem (cf. Meng et al., 2001; Liu and Wang, 2015) lies at that, in addition to physical links such as rail and road legs, which we refer to as carrier links throughout the rest of paper, the network planner needs to determine whether or not to improve/establish transfer links that represent cargo transfer processes between transportation modes/links at intermodal hubs. Thus, the network planner’s decisions are binary and four-fold: (i) establish a new carrier link (or not), (ii) expand the capacity of an existing carrier link (or not), (iii) establish a new transfer link, and (iv) expand the capacity of an existing transfer link. Note that locating a potential intermodal hub is determined at the same time with the decision on building new transfer links at that hub. If it is optimal to build at least one transfer link at a non-existing but candidate hub, the hub should be opened too; the hub location shall not be chosen otherwise. Since carriers and hub operators are final users of the network, the network planner is motivated to minimize the total operating cost of carriers and hub operators while making decisions on network design, in order to attract more business and retain more users. Meanwhile, there is a maximum budget that the network planner must not overspend when making decisions. The network planner would have to acquire the exact distribution of cargo flows over the network before she/he could move forth to calculate the total operating cost of carriers and hub operators under a proposal of network design. The network planner determines the network structure, but the cargo flow distribution is out of her/his control; instead, the cargo flow distribution is a result of the route choice of intermodal operators. Intermodal operators. Intermodal operators are decision makers who route cargo shipments through intermodal routes between each O–D pair using the services of carriers and hub operators along the routes, assigning cargo flows over the whole network. As a reward, the carriers and hub operators make revenue by earning service fees from intermodal operators. From the viewpoint of an analyst/modeler, the network planner could use a viable model to forecast the cargo flow distribution over the intermodal network. There exists a variety of models that have been used to predict the network flow distribution based on various assumptions on the choice behavior of route selectors. Classical studies on hub-and-spoke network design (Ernst and Krishnamoorthy, 1998; O’Kelly, 1987) assumed that unit transportation cost is constant on each link, no congestion effects (also known as the flow-time interaction that the unit travel/transfer time over each carrier/transfer link increases as more cargo flows are dealt with on the link) exist, and all cargo flows will be routed through the least-cost route between each O–D pair. Meng and Wang (2011) considered that congestion effects exist, each route’s utility (a weighted sum of transport time and cost, which measures the preferences of intermodal operators towards the route) is deterministic, and the route choice of intermodal operators follows the user-equilibrium principle that each intermodal operator would choose the largest-utility route, and at equilibrium, no intermodal operator could be better off by shifting to use another route. However, the assumption of deterministic utilities seems somewhat restrictive, especially when we consider a population of intermodal operators who have quite distinct perceptions on which route is of the largest utility. Discrete choice models, such as multinomial logit (MNL) model, nested logit model, and mixed logit model, can be used to fix the drawback (BenAkiva and Lerman, 1985; Train, 2003). Within the framework of discrete choice modeling, route utilities are considered to be random over the population of decision makers (i.e., intermodal operators in our study) and a decision maker chooses each route from a set of available routes with a certain probability. We will thus describe the choice behavior of intermodal operators with random route utilities using discrete choice models and the user-equilibrium principle that no intermodal operator would be better off by changing to choose another route at equilibrium. Unlike as is often assumed in the literature, we do not limit ourselves to a particular distribution for random route utilities and will use a fixed-point equation to model the cargo flow distribution over the network under a general discrete choice model. Such a treatment generalizes a broad family of cases where the choice behavior is described by a specific choice model such as the MNL model, nested logit model, or the mixed logit model.

78

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

1.1. Positioning and objective An intermodal network design problem with multiple stakeholders for cargo transportation was first introduced by Meng and Wang (2011). The authors formulated the problem as a mathematical program with equilibrium constraint (MPEC) model by assuming that all route utilities are deterministic. The MPEC model was then solved using a hybrid genetic algorithm, which could, however, possibly return a suboptimal solution. In addition, neither did the work offer a method to find a global optimum nor did it provide insights into how well the suboptimal solution performs compared with a globally optimal solution. Though it has been long since the work of Meng and Wang (2011), to the best of our knowledge, there exist no studies investigating the ITND problem under a general route choice model of intermodal operators with random route utilities. There exist no studies either that propose a solution algorithm to find globally optimal solutions to the ITND problem. Our study is aimed at filling the void and would make significant improvements over the work of Meng and Wang (2011). To be precise and concise, the objective of this work is to determine optimal decisions on whether or not to improve/build carrier and transfer links in an existing intermodal network to minimize the total operating cost of carriers and hub operators subject to the maximum budget expendable, accounting for the route choice behavior of intermodal operators with generally distributed random route utilities. 1.2. Challenges, our remedies, and contributions To achieve the objective, we formulate the ITND problem as a mixed-integer nonlinear programming model, called the ITND model, which incorporates three important features of the problem: (i) economies of scale for carriers and hub operators, (ii) congestion effects, and (iii) route choice behavior of intermodal operators. We use piecewise linear cost functions to capture the economies of scale, adopt flow-dependent time functions to reflect the congestion effects, and include a fixed-point equation to describe the route choice behavior of intermodal operators. Not surprisingly, the resulting ITND model is computationally challenging. The model is a mixed-integer nonlinear program with its objective function involving piecewise linear functions and with its constraints involving a fixed-point equation, making it very hard to solve. What makes it even harder is that the model will remain non-convex, even if we drop the integrality of all integer variables, due to the non-convexity of the fixed-point equation. To tackle the challenges, firstly, we re-formulate the ITND model as a variety of relaxed or equivalent but simpler mathematical programming models by exploiting the characteristics of piecewise linear cost functions and gradually relaxing constraints. We also provide theoretical results to justify the equivalences. Next, based on all these models, we solve the ITND problem using two strategies. For the first strategy, we are interested in developing an algorithm to find a globally optimal solution. For the second one, we will relax one of the equivalent models as a nonlinear optimization model and solve it using a traditional nonlinear optimization algorithm such as the steepest descent search method, trust-region algorithm, or interior-point method. This is the same as solving a continuous version of the ITND problem and generates a local optimal solution to the continuous ITND problem, in which, instead of determining whether or not to build/expand links, the network planner needs to determine how much capacity shall be added to each link. However, we might end up with getting fractional network design variables that cannot directly be applied to our (discrete) ITND problem. Thirdly, we propose two policies that transform possibly fractional solutions obtained from solving the relaxed nonlinear model to binary solutions. These solutions can be used as heuristic solutions to the hard ITND model. We also test the two solution strategies with numerical experiments with the hope to throw some light on the performance comparison between heuristic methods and a global optimization algorithm and on how to choose a proper solution approach for the ITND problem. Last but not least, while our model and algorithms are designed with the purpose of solving the discrete ITND problem, they could also be applied to solve the continuous ITND problem for both cargo and passenger transportation. Contributions. This work makes the following substantial contributions to the literature of intermodal transportation network design: (1) To the best of our knowledge, our work is the first to address the ITND problem that unifies three important features of multiple stakeholders: economies of scale, congestion effects, and a route choice model with generally distributed random route utilities. This lifts the assumption that route utilities are deterministic or have to follow a particular distribution. As a result of this remedy, our work can find interesting applications to a broad range of settings where the route behavior of intermodal operators is described by the MNL model, nested logit model, or mixed logit model. (2) Our work generalizes several assumptions that have been deemed appropriate in traditional hub-and-spoke network design problems (O’Kelly, 1987), but appears restrictive for intermodal cargo transportation: (i) our work involves cargo transfers between various modes of transport at intermodal hubs, (ii) it does not assume fully-connected hubs, and (iii) it does not require each route can at most traverse through two intermodal hubs. (3) We formulate the problem as a mixed-integer nonlinear programming model with a fixed-point constraint describing the route choice behavior of intermodal operators, i.e., the ITND model, and develop a variety of mathematical programs that are either equivalent to or relaxing the ITND model. We also provide theoretical results that validate the

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

79

equivalences between various models. In addition, we are able to show the existence and uniqueness of the solution to the fixed-point equation without restricting to strictly increasing flow-dependent time functions, which generalizes the results of traffic assignment models (Cascetta, 2009). (4) We present a nonlinear optimization algorithm to solve a relaxed nonlinear ITND model, from which we can obtain a local optimal solution to the continuous ITND problem as a byproduct. We introduce two policies that transform possibly fractional solutions obtained from solving the nonlinear model to integer solutions, which can be used as heuristic solutions to the ITND model. (5) We develop a branch-and-bound based global optimization algorithm to generate a globally optimal solution to the ITND model. The algorithm embeds a subroutine to obtain an upper bound and a subroutine to evaluate a lower bound by solving a mixed-integer linear program (MILP) obtained from relaxing the ITND model. (6) We test the heuristic solutions and the global optimal solution with numerical experiments and shed some light on how well heuristic solutions perform compared with the global optimizer. 2. Literature review We classify the literature related to our research into three branches: (1) hub-and-spoke network design, (2) network design for intermodal cargo transportation, (3) network design under route choice models. 2.1. Hub-and-spoke network design Crainic and Kim (2007) commented that an intermodal transportation network can be naturally organized as a hub-andspoke network. The hub-and-spoke network design problem focuses on locating a certain number of hubs from a set of potential locations. As one of the earliest works on hub-and-spoke network design, O’Kelly (1987) solved a p-hub median problem and applied the results to passenger airline network design. A great deal of literature showed up afterward in the literature. Interested readers might refer to Alumur and Kara (2008) for quite a complete and inspiring literature review of the topic. The studies on hub-and-spoke network design often share several common assumptions (Ernst and Krishnamoorthy, 1998), for example, hub link cost is discounted by a factor less than one to reflect the economies of scale, each route can at most pass two hubs, and hubs are fully interconnected. These assumptions are made in part due to the practical characteristics of passenger airline networks, for which the hub-and-spoke network design problem finds the very first applications, and partly for easing the model formulation or solution design. Luna and Mahey (20 0 0) considered a global optimization of capacity expansion and flow assignment problem and formulated the problem as a non-convex mathematical program. The authors solved the hard non-convex problem by convexifying arc functions and showed how to calculate the bounds for the approximation. However, it becomes more complicated when we talk about intermodal transportation and the assumptions above will need to be relaxed. Some of these assumptions have been relaxed. O’Kelly and Bryan (1998) and Kimms (2006) considered a more general flow-dependent cost function to describe the economies of scale over hub links. Elhedhli and Hu (2005) incorporated congestion effects in the hub-and-spoke network design. Our work generalized the assumptions above, and furthermore, we also add the features of multiple stakeholders and the behavior of route choice makers, which have not thus far appeared in the studies on the classical hub-and-spoke network design problem or its variants. 2.2. Network design for intermodal cargo transportation There exists an emerging trend of studies on intermodal network design for cargo transportation. Alumur et al. (2012) formulated the multimodal hub location problem as an MILP, in which constraints on delivery times are enforced and a variety of decisions are made, including hub locations, spoke-hub allocations and hub link modes. In a subsequent paper, Serper and Alumur (2016) added a new feature to determine the number of vehicles of each type to be deployed on the hub network in addition to hub locations and spoke-hub allocations. Demir et al. (2016) studied an energy-ware intermodal service network design problem with uncertainties in travel time and demands. Meraklı and Yaman (2016) proposed a robust mixed-integer linear programming model for the p-median hub location problem with a polyhedral uncertainty set of demands. These studies, however, considered a single stakeholder and involved no route choice behavior or congestion effects. 2.3. Network design under route choice models There is another stream of literature on network design with the choice behavior of route selectors, with major applications to urban passenger transportation rather than cargo transportation. LeBlanc (1975) applied a branch-and-bound algorithm to solve a bi-level programming model for urban transportation network design, in which the lower-bound problem was solved by assuming deterministic route utilities and the userequilibrium principle. Meng et al. (2001) considered a continuous network capacity design problem for urban transportation with the user-equilibrium principle and solved the problem using an augmented Lagrangian algorithm. Yamada et al. (2009) established a bi-level programming model for freight transportation network design, in which the lower-level problem is a multimodal multi-class user-equilibrium traffic assignment problem with deterministic route utilities. A series of

80

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

heuristic approaches including a genetic algorithm and tabu search based procedures were tested in realistic networks. As discussed, Meng and Wang (2011) introduced an intermodal network design problem for cargo transportation with deterministic route utilities and solved the problem suing a hybrid genetic algorithm. Wang et al. (2013) proposed a global optimization algorithm to solve a discrete transportation network design problem using the relation between the userequilibrium and system optimum models for traffic assignment. These studies, however, did not consider a full set of features of multiple stakeholders, economies of scale, and random route utilities. Loureiro and Balston (1996) looked into a bi-level multi-commodity network design problem for determining investment policies for intercity freight transportation networks by assuming that the behavior of intermodal operators in route choice follows the MNL model. A two-step heuristic algorithm based on column generation was designed to solve the model. In a recent work of Liu and Wang (2015), a global optimization algorithm was proposed to solve a continuous network design problem with the MNL describing the behavior of choice makers under Gumbel-distributed route utilities. These studies did not consider multiple stakeholders or economies of scale. In addition, the widely chosen MNL model bears the well-known independence from irrelevant alternatives (IIA) (Luce, 1959, p. 9) property that the ratio of choice probabilities of any two alternatives is independent of the presence or absence of a third alternative in the choice set. The IIA property, on one hand, could be exploited to ease the problem; on the other hand, cannot model the preference heterogeneity among choice makers, for which the nested logit model or mixed logit model seems a better choice. As mentioned, our work does not assume a particular route choice model and can be applied with a family of choice models, including the MNL, nested logit model, and mixed logit model. As an example, we use the nested logit model to describe the route choice behavior in our numerical tests. With all the explanations above, we shall devote Section 3 to representing an intermodal network as an operational network with transfer links. We then pass on to Section 4 to introduce various stakeholders involved and define the ITND problem. In the succeeding Section 5, we formulate the problem as a mixed-integer nonlinear programming model, i.e., the ITND model, and propose a variety of mathematical programs either equivalent to or relaxing the ITND model. Next, we present in Section 6 a nonlinear optimization algorithm and a global optimization algorithm to solve the problem. A very important Section 7 will then be given to test the two solution algorithms with numerical experiments. In the final Section 8, we provide a brief recapitulation of the whole paper and offer some concluding remarks. The technical proofs, supplemental materials, and test data are provided in Appendix A–Appendix C. 3. Intermodal network model We introduce in this section the representation of the physical intermodal network and discuss how to transform it into an operational network, for which the planner makes network design decisions. The network representation is similar to but different from that described in Meng and Wang (2011) in two aspects: (i) We involve no-choice links to reflect the operational reality that there exist intermodal operators who possibly choose not to transport cargo through the network, while Meng and Wang (2011) did not. This adds another dimension of flexibility to the model and is more realistic. (ii) We consider the feasible sets of route and link flows to be independent of network design decisions. In Meng and Wang (2011), both the sets are defined as a function of network design decisions, which complicates the model by incorporating the interactions between network flows and network design decisions into the flow sets. However, with the network representation in this work, the sets are defined to be independent of network design decisions and we leave the interactions to be described by the fixed-point equation that models the choice of intermodal operators. By doing so, we obtain a simpler model with a more compact structure. Before getting into the formal problem statement, we first state the notational conventions in this paper. We consider all the vectors as column vectors unless otherwise stated. However, for the economical use of space, we use x = (x1 , x2 , . . . , xd ) inside the text to describe a vector even though x ∈ Rd is a column vector, where d ∈ N. Let x1 and x2 denote the 1-norm and 2-norm of any vector x ∈ Rd , where d ∈ N. Let |A| denote the cardinality of any set A. Define x+ := max{x, 0} for all x ∈ R. We interpret ∞ to be an ideally large number with the following arithmetic operations: 1/∞ = 0, 1/0 = ∞, e−∞ = 0, −∞ + a = −∞ for all a ∈ R, and ∞ − ∞ = 0. By X ∼ F, we mean that any random variable X follows distribution F. Physical network. Consider a physical intermodal network under consideration of re-design. Let N denote the set of nodes in the network, including spoke nodes (without cargo transferring), existing intermodal hubs, and potential intermodal hubs to be established. Denote by M the set of modes. For example, M = {rail, truck, ocean}. There may exist a link (or arc, leg) pointing from a node i ∈ N to node j ∈ N with mode m ∈ M, where i is the tail and j is the head of the link. Note that we allow the existence of multiple parallel modes between i and j and each link is specified by its head, tail, and mode. Let Ac ⊂ {(i, j, m ) : i, j ∈ N, m ∈ M} denote the set of existing and potential bi-directional links between nodes. On each link, there exists a carrier that provides cargo transportation services utilizing vehicles, crew, and equipment. This is why we refer to these links as carrier links and Ac as the set of carrier links. For any node n, let Ain (n ) ⊂ Ac denote the set of carrier links pointing into node n and Aout (n ) ⊂ Ac be the set of carrier links stemming out of node n. Thus, the network can be represented by a directed graph G := (N, Ac ). Intermodal hubs and transfer links. Let H ⊂ N denote the set of existing and potential intermodal hubs, where containerized cargo is processed to be transferred from one mode (carrier link) to another mode (carrier link). For each intermodal hub h ∈ H, we represent the transfer process from a link in Ain (h ) to another link in Aout (h ) as a transfer link at the

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

81

Fig. 1. Intermodal network representation: physical and operational networks.

hub. The operational process over a transfer link is implemented by the intermodal hub operator using crew and equipment such as cranes and forklifts. We now describe how to generate transfer links at each intermodal hub h ∈ H. Note that all carrier links in Ain (h ) share the same head, i.e., hub h, and all carrier links in Aout (h ) originate from the same tail, i.e., hub h. We make |Ain (h )| + |Aout (h )| copies of hub h and replace the hub by these duplicate nodes. For each carrier link a ∈ Ain (h ), we assign a copy of hub h as its head. For each carrier link a ∈ Aout (h ), we allocate a copy of hub h to be its tail. The assignment is executed in a way that each carrier link in Ain (h ) ∪ Aout (h ) is assigned a unique copy hub and each copy is allocated to a unique carrier link. Let Nin (h ) denote the set of heads (copies) of all carrier links in Ain (h ) and let Nout (h ) denote the set of tails (copies) of all carrier links in Aout (h ) after assignment. Thus, the set of transfer links at hub h is denoted by A(h ) := {(i, j ) : i ∈ Nin (h ), j ∈ Nout (h )} and the set of all transfer links is represented by At := ∪h∈H A(h ). Operational network. Let N := ∪h∈H (Nin (h ) ∪ Nout (h ) ) ∪ N denote the set of nodes after transfer link generations, which is obtained by augmenting N with the copies of intermodal hubs. Intermodal operators send containerized cargo from origin nodes to destination nodes. Let O ⊂ N be the set of origin nodes and D ⊂ N be the set of destination nodes. Denote by W ⊂ O × D the set of origin-destination (O–D) pairs. We associate each O–D pair w ∈ W with a no-choice link 0w , which represents the action of intermodal operators not choosing any route between w and is considered as a direct dummy link connecting w. Let Ad := {0w , w ∈ W } denote the set of no-choice dummy links. Define A := Ac ∪ At ∪ Ad and we denote by G := (N , A ) the operational network generated from G = (N, Ac ). Note that there exists at least one route between each O–D pair. Fig. 1 shows the example of a physical intermodal network in Fig. 1a transformed into an operational network in Fig. 1b, in which intermodal hub h = 3 in G is replaced by its duplicate nodes, 5, 6, 7, and 8 in G. The figure also shows explicit values for above-defined various sets. Feasible sets of route and link flows. For any w ∈ W, denote by Rw the set of intermodal routes connecting w (including a,r 0w ∈ Rw ) and let δw be one if route r ∈ Rw traverses a ∈ A; zero otherwise. Note that the no-choice link 0w exists as both 0w ,r a link and an intermodal route connecting w. Link 0w cannot be part of any route r = 0w . By notation, δw = 1 if r = 0w ; a,r it is zero otherwise. Denote by w := (δw , a ∈ A, r ∈ Rw ) the link-route incidence matrix for each w ∈ W, where links and routes correspond to rows and columns of w , respectively. Let R := ∪w∈W Rw be the set of all routes between all O–D pairs. Denote by  := (w , w ∈ W ) the link-route incidence matrix for graph G, where columns are ordered by O–D pairs. Let Dw ∈ R+ := [0, ∞ ) denote the demand that needs to be transported between O–D pair w ∈ W. For any link a ∈ A, denote by va ∈ R+ the cargo flow representing number of standardized cargo units (TEUs or FEUs) loaded on a. Let v := (va , a ∈ A ) be the vector of link flows. For each w ∈ W and r ∈ Rw , let fw,r ∈ R+ denote the cargo flow loaded on route r. Define fw := ( fw,r , r ∈ Rw ) as the vector of route flows for all w ∈ W and denote by f := (fw , w ∈ W ) the vector of route flows over network G. Let v be the feasible set of link flows defined as

  v := v ∈ R|A| : ∃ f ∈ R|+R| s.t. v = f, fw 1 = Dw , ∀ w ∈ W .

Let f be the feasible set of route flows, which is defined as,

   f := f ∈ R|+R| : fw 1 = Dw , ∀ w ∈ W .

82

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Proposition 1 shows that both the sets of link and route flows are nonempty. Proposition 1. It holds that v = ∅ and  f = ∅ and both v and f are compact sets. 4. Problem definition We discuss in this section the roles of various stakeholders, including the network planner, carriers, hub operators, and intermodal operators, and formally define the ITND problem. 4.1. Decision variables of the network planner Denote by A0 ⊂ A the set of existing links that do not require any actions such as maintenance or rehabilitations. Clearly, it follows that the set of no-choice links Ad ⊂ A0 . The network planner may need to expand the capacities of some of the existing links or establish new links to eliminate traffic congestions. Let A1 ⊂ A denote the set of existing links to be improved. Let A2 ⊂ A denote the set of potential links to be established. We have that Ai ∩ A j = ∅ for all i, j ∈ {0, 1, 2}, i = j and A = ∪2i=0 Ai . Note that if A0 = A1 = ∅, the network planner is faced a problem to establish a brand new network. Define the decision variables of the network planner as follows:



xa =

1 0

if link a is established or its capacity is expanded, otherwise,

∀ a ∈ A.

Thus, the set of decision variables is defined as,



X := x ∈ {0, 1}|A| : xa = 1,

 ∀ a ∈ A0 .

Any vector x ∈ X specifies a feasible design of the network’s structure and is referred to as a network design vector. Define conv(X ) := x ∈ [0, 1]|A| : xa = 1, ∀ a ∈ A0 as the convex relaxation of X . Let a ∈ R+ denote the original capacity of link a ∈ A, where  a > 0 for a ∈ A0 ∪ A1 and a = 0 for a ∈ A2 . Let γa ∈ R+ denote the expansion of the capacity of link a ∈ A such that γ a > 0 for a ∈ A1 ∪ A2 and γa = 0 for a ∈ A0 . Thus, the link capacity is a (xa ) := a + xa γa , i.e., if the link is established or its capacity is expanded (xa = 1), the link capacity will be updated up to a + γa ; the capacity remains at  a otherwise (xa = 0). Let Ba ∈ R+ be the budget invested in expanding or establishing link a ∈ A with Ba = 0 for all a ∈ A0 . Let H0 ⊂ H denote the set of existing intermodal hubs and H1 := H \ H0 denote the set of potential intermodal hubs. It follows that ∪a∈A(h ) ⊂ A2 if h ∈ H1 is a potential hub. Note that if h ∈ H1 is a potential hub, we consider that the total cost for building h can be approximated by summing up the costs for building transfer links within that hub. Under such a consideration, only those transfer links that are determined to be built (by the optimal network design vector x) will be counted toward the total cost. The transfer links not chosen to be established will not contribute to the total cost. Thus,  a∈A (h ) xa Ba will be resulting total cost for building the hub. Remark 1. It is worthwhile to note that locating a potential hub is determined by the optimal decisions on building up the  transfer links at the hub, that is to say, if x is the optimal network design vector and a∈A(h ) xa ≥ 1, then h is opened; it is not chosen to be a hub otherwise. Remark 2. It holds that X = ∅, since x ∈ R|A| with xa = 1 for all a ∈ A resides in X . Then, X ⊂ conv(X ) = ∅. Remark 3. It follows from Proposition 1 and Remark 2 that v × conv(X ) = ∅. 4.2. Intermodal operators’ route choice model 4.2.1. Intermodal route utilities Transport fare. In general, the service provider over a carrier or transfer link charges its customers, i.e., intermodal operators, a predetermined fare for transporting/transferring one unit of cargo, e.g., one TEU, through the link (CRRC, 2012). We use pa ∈ R+ to denote this fare amount and consider it to be a constant for any a ∈ Ac ∪ At . Travel/transfer time. Denote by ta the time function for transporting/transferring one unit of cargo, e.g., one TEU, on link a ∈ Ac ∪ At . We consider flow-dependent unit cargo handling time over a link to capture the congestion effect for any a ∈ Ac ∪ At . More specifically, let ta : R+ × [0, 1] → R+ ∪ {∞} be the extended real-valued time function that represents the travel/transfer time for per unit of cargo over link a ∈ Ac ∪ At . We make the following assumption for time functions. Assumption 1. For each a ∈ Ac ∪ At and any xa ∈ [0, 1], we assume that ta ( ·, xa ) is a continuous and monotone nondecreasing function. Note that, not as is often assumed in the literature (Cascetta, 2009, pp. 311–312), we do not require ta (va , xa ) to be strictly increasing in va for any xa ∈ [0, 1].

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

83

For each carrier link a ∈ Ac , we consider the following travel time,

ta (va , xa ) := ta0 +

θa0 ∀ a ∈ Ac , (a (xa ) − va )+

(1)

with the understanding that 1/0 = ∞, where ta0 ∈ R+ is a constant representing the congestion-free travel time once the link is built/improved, the second term captures the congestion effect over the link, and θa0 ∈ R+ is a constant parameter to be estimated from data. For each transfer link a ∈ At , we use the following cargo transfer time:

ta (va , xa ) :=

θa0 ∀ a ∈ At , (a (xa ) − va )+

(2)

subject to the convention that 1/0 = ∞, where θa0 ∈ R+ is a constant parameter to be estimated from data. Functions (1) and (2) are formulated based on Kingman’s approximation of the cargo processing time under congestion, where the cargo handling process is thought to be a queueing system, va is interpreted as the arrival rate of container flows, and  a (xa ) is explained as the service rate of the link. Note that, we allow the cargo travel/transfer time over a link to be infinite, when a link ends up being unestablished or completely congested due to the traffic density va / a (xa ) exceeding unity. ¯ := R ∪ {−∞} as its systematic utility function, which Systematic utility. For each a ∈ Ac ∪ At , define ua : R+ × [0, 1] → R is represented by

ua (va , xa ) := −α pa − β ta (va , xa ),

(3)

where α , β ∈ (0, ∞) represent price and time sensitivity of intermodal operators, both α and β are considered to be identical across the entire population of intermodal operators, and need to be estimated from data. For any O–D pair w ∈ W and for each intermodal route r ∈ Rw \ {0w }, define gw, r as the systematic utility of r representing the preference of intermodal operators towards intermodal route r, which is of the additive form,

gw,r (v, x ) :=



δwa,r ua (va , xa ), ∀ r ∈ Rw , w ∈ W.

(4)

a∈A

For no-choice link (route) 0w , let u0w = gw,0w be a constant representing its systematic utility. Note that gw,0w is free of decision variable x and cargo flow v. In practical applications, an intermodal operator would give up transporting cargo if the utilities of all the real available routes are lower than a sufficiently small level. Thus, it is reasonable to assume |gw,0w | < ∞, which is equivalent to that the no-choice route has a nonzero chance to be chosen for each O–D pair. Define the vector of link systematic utilities as u(v, x ) := (ua (va , xa ) : a ∈ A ) for all (v, x ) ∈ v × conv(X ). The above relation (4) indicates that the systematic utility of a route is the summation of systematic utilities of the links that the route traverses through. For each w ∈ W, define gw (v, x ) := (gw,r (v, x ), r ∈ Rw ), where it holds that,

g w ( v , x ) =  w u ( v, x ) ,

∀ (v, x ) ∈ v × conv(X ).

4.2.2. Random utility based choice model Random utility. An intermodal operator desires her cargo to be transported from the origin to destination of an O-D pair w ∈ W and needs to make a selection of routes from Rw to accomplish the cargo transportation. As discussed in Section 1, we consider randomly distributed route utilities. Let ξ w, r be the continuously distributed random error added to the systematic utility gw, r with mean zero and support w,r ⊂ R. Define ξ := (ξw,r , w ∈ W, r ∈ Rw ). Let  be the support of ξ . It follows that  = w∈W r∈Rw wr . Let P denote the probability measure defined on . For any w ∈ W, the random route utility g˜w,r of route r ∈ Rw is defined as

g˜w,r (v, x, ξ ) := gw,r (v, x ) + ξw,r ,

∀ w ∈ W, r ∈ Rw .

Choice operator. Given (v, x ) ∈ v × conv(X ), the probability of intermodal operators choosing r ∈ Rw , where w ∈ W, is ¯ |Rw | → [0, 1], which is defined by given by the choice function qrw : R



qrw



(gw (v, x )) := P r ∈ arg max g˜w, j (v, x, ξ ) j∈Rw

 

.

(5)

¯ |Rw | for all w ∈ W and r ∈ Rw . Assumption 2. We assume that qrw is a differentiable function on domain R This assumption holds under a variety of cases depending on the distribution of ξ . For example, if (ξw,r , r ∈ Rw ) is a vector of independent Gumbel random variables, qrw is of a closed-form expression, known as the MNL model, and it is continuously differentiable. Other examples include nested logit model, probit model, and mixed logit model (Cascetta, 2009). ¯ |Rw | → [0, 1]|Rw | by qw (gw (v, x )) := (qrw (gw (v, x )), r ∈ Rw ) for all w ∈ W. Define the choice operator qw : R

84

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Remark 4. The probability qrw (gw (v, x )) in (5) is well defined in the sense that qw (gw (v, x ))1 = 1, since it follows from the additive rule of probability measure P that,

qw (gw (v, x ))1 =





j∈Rw

r∈Rw



=P



P r ∈ arg max{g˜w, j (v, x, ξ )}

r∈Rw





r ∈ arg max{g˜w, j (v, x, ξ )} j∈Rw

= 1.

Remark 5. It follows from |gw,0w | < ∞ that qrw (gw (v, x )) = 0 if gw,r ((v, x )) = −∞ for any r ∈ Rw and w ∈ W. Cargo flow distribution and fixed-point equations. According to the user equilibrium principle with random route utilities (Cascetta, 2009), the cargo flow distributed over routes is given by the following route-version fixed-point equation

fw = Dw qw (gw (f, x )),

∀ w ∈ W, f ∈  f , x ∈ conv(X ),

(6)

and that distributed over links is given by the link-version fixed-point equation,

v = F (v, x ) :=



Dw w qw (gw (v, x )),

(7)

w∈W

where F : v × conv(X ) → v is called the fixed-point operator. 4.2.3. Characterization of the solution to the link-version fixed-point equation Proposition 2 states that the fixed-point operator F maps any (v, x) to set v . Proposition 2. For any (v, x ) ∈ R|A| × conv(X ), it holds that F(v, x) ∈ v . Lemma 1 shows the monotonicity of the choice operator qw . Not as is conventionally considered, qw is defined on an ¯ |Rw | := r∈R [R ∪ {−∞}] in our paper. We show that the monotonicity of qw still exists on the extended extended domain R w domain subject to our notational conventions. Lemma 1. For each w ∈ W, the choice operator qw is a monotone operator, i.e.,

(qw (g1 ) − qw (g2 )) (g1 − g2 ) ≥ 0, ∀ g1 , g2 ∈ R¯ |Rw | . We now discuss the existence and uniqueness of the solution to the fixed-point equation (7). Note that, under Assumption 1, we drop the assumption that time function ta ( ·, xa ) has to be strictly increasing, which has long been petted in network flow assignment (Cascetta, 2009, pp. 311–312), partly for simplifying the problem, but does not hold in our study. Without confining time functions to strictly increasing maps, however, we are still able to show the uniqueness of the solution to the fixed-point equation (7) in Proposition 3, and this pushes forth the state of the art of traffic assignment theories. Proposition 3. For any x ∈ conv(X ), there exists a unique v ∈ v such that v = F (v, x ). Proposition 4 states that if a new link a ∈ A2 is not chosen to be established, then there exist no flows loaded over it. Meanwhile, all intermodal routes that traverse through a have zero flows. Proposition 4. Consider (v, x ) ∈ v × conv(X ) that solves the fixed-point equation (7). If xa = 0 for any a ∈ A2 , i.e., link a is not established, then the following holds: a,r (1) fw,r = 0 for all r ∈ Rw and w ∈ W such that δw = 1, (2) va = 0.

4.3. Economies of scale and cost functions for carriers and hub operators Cargo transportation/transfer operations incur costs to carriers/hub operators, which could include installation cost, vehicle operating cost, crew salaries, and management cost. Let Ca : R+ → R+ denote the cost function for the carrier/hub operator that operates link a ∈ Ac ∪ At . Therefore, Ca (va ) is the cost incurred for handling cargo flow of va . It has been widely observed that the marginal cost decreases with the increasing amount of cargo that a carrier/hub operator handles. This size advantage leads to the so-called economies of scale. The cost function Ca is often assumed to be a nondecreasing concave function to reflect the economies of scale. We consider that Ca (va ) is a piecewise linear function of a minimization form,

Ca (va ) := min{bia va + eia }, i∈Ia

∀ a ∈ Ac ∪ At ,

where Ia := {1, 2, . . . , Ia } and Ia ∈ N is the number of linear segments.

(8)

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

85

Fig. 2. The piecewise linear cost function.

Fig. 2 shows the piecewise linear cost function Ca for link a, in which b1a , b2a , . . . , bIaa are slopes, e1a , e2a , . . . , eIaa are the corresponding intercepts, and a0 , a1 , a2 , . . . , aIa −1 are the dividing points for all linear segments along the horizontal axis of flow va , where a0 = 0. Define aIa := ∞. Property 1. The following properties hold for cost function (8): b1a > b2a > . . . > bIaa ≥ 0, representing the decreasing marginal cost and reflecting the economies of scale. 0 ≤ e1a < e2a < . . . < eIaa . Ca (0 ) = e1a ≥ 0. For any a ∈ Ac ∪ At and va ∈ R+ , the minimum in (8) is attained at only one segment i ∈ Ia , if va ∈ (ai−1 , ai ), or two segments i − 1 and i if va = ai−1 . (5) Ca (va ) is a concave function. (1) (2) (3) (4)

Note from Property 1(3) that we allow a fixed installation cost associated with each link, which represents the possible capital cost that a carrier/hub operator invests even though there are no cargo flows handled for business. For each a ∈ Ac ∪ At and va ∈ R+ , denote by ia ∈ Ia such that va ∈ [aia −1 , aia ) throughout the remaining paper. We now present Lemma 2 that provides fundamental but useful results for any real-valued concave function of a single variable. Lemma 2. (Subadditivity) Let f : R → R be a concave function that is not a zero function. Then the following holds: (1) If f(0) > 0, we have yf(x) > xf(y) for any x, y ∈ R+ such that x < y. (2) If f(0) > 0, we have f (x ) + f (y ) > f (x + y ) for any x, y ∈ R+ . (3) If f (0 ) = 0, we have yf(x) ≥ f(y)x for any x, y ∈ R+ such that x < y, where the inequality holds as equality if and only if x = 0 or f is a linear function passing through 0. (4) If f (0 ) = 0, we have f (x ) + f (y ) ≥ f (x + y ) for any x, y ∈ R+ , where the inequality holds as equality if and only if xy = 0 or f is a linear function passing through 0. Lemma 2 and Property 1(3) imply that cost function Ca (va ) is subadditive in flow variable va . Corollary 1. For any x, y ∈ R and all a ∈ Ac ∪ At , if e1a > 0, then Ca (x ) + Ca (y ) > Ca (x + y ) and if e1a = 0, then Ca (x ) + Ca (y ) ≥ Ca (x + y ). We now formally state the ITND problem: the network planner needs to determine which links to be established/improved to minimize the total operating cost for carriers and hub operators, subject to the network flow distribution resulting from intermodal operators’ route choice. 5. Model Below, we formulate the ITND problem as a mixed-integer nonlinear programming model, which we call the ITND model. The ITND model involves both piecewise linear costs function to capture the economies of scale and a fixed-point constraint to reflect intermodal operators’ route choice behavior, making it difficult to solve. We will, therefore, develop a variety of mathematical programming models that are either equivalent to or relaxing the ITND model, in hopes of reducing the hardness of the problem.

86

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

5.1. The ITND model We formulate the ITND problem with intermodal operators’ route choice behavior as the following ITND model,

(ITND) min C (v, x ) := v,x

a∈Ac ∪At



subject to:



Ca (va ) +



(xa − 1 )e1a + ρ

a∈A2



va

(9)

a∈Ad

Ba xa ≤ B,

(10)

a∈A

va ≤ a (xa ) ∀ a ∈ Ac ∪ At ,

(11)

v = F ( v, x ) ,

(12)

( v , x ) ∈ v × X ,

(13)

where ρ ∈ R+ is a positive number that represents the penalty on losing cargo flows to no-choice links, which the network planner tries to hedge against. The minimization of objective function (9) reflects the network planner’s will to minimize the sum of the network construction cost and the total operating cost for carriers and hub operators over network G. The network planner is motivated to do so to retain more users and businesses. Inequality (10) is the budget constraint. We use Constraints (11) to enforce that a link cannot process cargo flows more than its capacity. Fixed-point equation (12) is included as a constraint for two purposes: (i) it models intermodal operators’ route choice behavior and distributes cargo flows over the intermodal network for a given candidate network design vector x ∈ X and (ii) it preserves the network flow balance relations between links and intermodal routes. Constraint (13) simply specifies the feasible sets for v and x.  Remark 6 is obtained based on the inclusion of the second term “ a∈A2 (xa − 1 )e1a ” into the objective value (9). Remark 6. Consider any feasible solution (v, x) to the ITND model and consider any a ∈ A2 to be established. If xa = 0, i.e., link a is not established, then the operating cost Ca (va ) for link a, including the fixed installation cost, is not counted towards the objective value. If xa = 1, then the operating cost, including the fixed installation cost, is counted into the objective value. To see this, if xa = 0 for any a ∈ A2 , it follows from Proposition 4(2) that va = 0. Then, it follows from Property 1(3) that the operating cost for link a is computed as Ca (0 ) = e1a , which cancels out with (xa − 1 )e1a = −e1a , in the second term of the objective value. If xa = 1, then the second term (xa − 1 )e1a = 0. This shows the remark. Remark 7 asserts that the objective function is continuous and it attains the minimum over the feasible set to the ITND model. Remark 7. The objective function C(v, x) is continuous and concave, since Ca (va ) is a concave continuous function. Since X is finite-dimensional, it follows from Proposition 3 that the feasible set to the ITND is a discrete compact set. Thus, continuous function C(v, x) attains its minimum over the compact feasible set. Proposition 5 asserts that capacity constraint (11) is redundant. Proposition 5. If (v, x) satisfies the fixed-point constraint (12), then it satisfies (11). 5.2. An equivalent mixed-integer nonlinear programming model The objective function of the ITND model involves piecewise linear cost functions, which brings up challenges for solving the model using a general optimization solver or designed algorithm. Inspired by the work of Croxton et al. (2003), we ease the problem by transforming the objective function into a linear function with complementary variables. For each a ∈ Ac ∪ At and i ∈ Ia , define yia ∈ {0, 1} as a segment variable such that



yia =

1 0

if va ∈ [ai−1 , ai ), otherwise,

and define via ∈ R+ as the flow variable such that



i a

v =

va 0

if va ∈ [ai−1 , ai ), otherwise.

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104







87



Define ya := yia , i ∈ Ia for all a ∈ Ac ∪ At , y := (ya , a ∈ Ac ∪ At ), and define z := v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At . Define the following cost function,

C1 (z ) :=

Ia   



bia via + yia eia +

a∈Ac ∪At i=1



(xa − 1 )e1a + ρ

a∈A2



va .

(14)

a∈Ad

We formulate the following ITND-1 model:

(ITND-1) C1∗

:=

z

subject to:

Ia

i=1

(10 ), (12 ), Ia i v = va i=1 a

∀ a ∈ Ac ∪ At ,

(15)

∀ a ∈ Ac ∪ At ,

yia = 1

ai−1 yia ≤ via ≤ ai yia yia ∈ {0, 1}

C1 (z )

min

(16)

∀ i ∈ {1, 2, . . . , Ia }, a ∈ Ac ∪ At ,

(17)

∀ i ∈ {1, 2, . . . , Ia }, a ∈ Ac ∪ At ,

(18)

via ≥ 0, ∀ i ∈ {1, 2, . . . , Ia }, a ∈ Ac ∪ At

(19)

via ≤ yia D ∀ i ∈ {1, 2, . . . , Ia }, a ∈ Ac ∪ At ,

(20)

x ∈ X,

(21)



where D := w∈W Dw and Constraints (15)–(19) exist to model the piecewise linear cost function for each link a. Constraint (20) is redundant in the ITND-1 model, but adding it will help reduce the problem size, as later showed by Theorem 2. Theorem 1 justifies the equivalence between the ITND model and ITND-1 model. Theorem 1. If z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) is an optimal solution to the ITND-1 model, then, (v, x) is an optimal solution to the ITND model. 5.3. An equivalent model that drops Constraints (16) and (17) Inequalities (15)–(19) involve a large number of constraints and integer variables (i.e., y), which seems to be complicating the ITND-1 model. Thus, we are interested in looking for possible ways to simplify the ITND-1 model. Define the following cost function,

C˜1 (z ) :=





e1a +

a∈Ac ∪At

Ia  

bia via + yia (eia − e1a )

 

+



(xa − 1 )e1a + ρ

a∈A2

i=1



va .

a∈Ad

We obtain an ITND-1R model by removing Constraints (16) and (17) from the ITND-1 model, which is presented below:

(ITND-1R ) C˜1∗ := min C˜1 (z ) z

subject to: (10 ), (12 ), (15 ), (18 ) − (21 ). For each a ∈ Ac ∪ At , define the adjusted cost function





C˜a (va ) := min bia va + eia − e1a , i∈Ia

∀ a ∈ Ac ∪ At .

Remark 8. It follows that C˜a (0 ) = 0 and that all properties in Property 1 are valid for function C˜a with eia replaced by eia − e1a for all i ∈ Ia . Theorem 2 asserts that Constraints (16) and (17) can be dropped from the ITND-1 model, which largely cuts down the problem size.

88

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Theorem 2. Consider any z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) feasible to the ITND-1R model and create zˆ = (v, x, yˆ , vˆ ia , i ∈ Ia , a ∈ Ac ∪ At ), where



yˆia

=



if va ∈ [ai−1 , ai ), otherwise,

1 0

ˆ ia

v =

va 0

if va ∈ [ai−1 , ai ), otherwise,

for each a ∈ Ac ∪ At . The following holds: (1) (2) (3)

zˆ is feasible to the ITND-1R model and C˜1 (z ) ≥ C˜1 (zˆ ). If z is optimal to the ITND-1R model, then zˆ is also optimal to the ITND-1R model and (v, x) is optimal to the ITND model. C˜1∗ = C1∗ .

Theorem 2 (1) suggests an intriguing result that, given (v, x), any feasible solution z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) is dominated by a feasible solution created as zˆ = (v, x, yˆ , vˆ ia , i ∈ Ia , a ∈ Ac ∪ At ). 5.4. An equivalent model relaxing integer segment variables Define the following cost function, Ia  

C2 (z ) :=





yia bia va + eia +

a∈Ac ∪At i=1



(xa − 1 )e1a + ρ

a∈A2



va .

(22)

a∈Ad

We also provide the following mixed-integer nonlinear program, called the ITND-2 model, which further relaxes the integrality of y,

(ITND-2)

C2∗ := min C2 (z ) z

∀ i ∈ {1, 2, . . . , Ia }, a ∈ Ac ∪ At , (10 ), (12 ), (16 ), (21 ).

subject to: yia ≥ 0

(23)

The non-convexity of the ITND-2 model stems from the route choice model (12) and the first term of the objective value. Thus, there may exist multiple local optimal solutions to the ITND-2 model. Proposition 6 characterizes any (local or global) optimal solution to the ITND-2 model. Recall that ia ∈ {1, 2, . . . , Ia } is such that va ∈ [aia −1 , aia ) for each a ∈ Ac ∪ At . Proposition 6. If z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) is a local (global) optimal solution to the ITND-2 model, then for each a ∈ Ac ∪ At , the following two cases can happen: (1) If va ∈ (aia −1 , aia ), then yiaa = 1 and yia = 0 for all i = ia . (2) If va = aia −1 , then yia = 0 for all i ∈ / {ia − 1, ia }. Proposition 6 and Property 1(4) imply Corollary 2. Corollary 2. If z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) is a local (global) optimal solution to the ITND-2 model, then  Ia i  i y b v + eia = biaa va + eiaa for all a ∈ Ac ∪ At , where, as defined, ia ∈ Ia such that va ∈ [aia −1 , aia ). i=1 a a a Corollary 2 yields Theorem 3 that justifies the equivalence between the ITND-2 model and ITND model. Theorem 3. If z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) is an optimal solution to the ITND-2 model, then, (v, x) is optimal to the ITND model and C2∗ = C1∗ . 5.5. An MILP model that relaxes the ITND model We formulate below an MILP model that relaxes the ITND-2 model and provide Theorem 4 to justify the relaxation.

(ITND-2R)

C˜2∗ := min C˜1 (z ) z

subject to: (10 ), (11 ), (15 ), (19 ) − (21 ), (23 ) v ∈ v .

(24)

Theorem 4 implies that the optimal objective value of the ITND-2R model is a lower bound on the optimal objective value of the ITND-2 model and the optimal objective value of the ITND-1R model. Theorem 4. It holds that C˜2∗ ≤ C2∗ = C˜1∗ .

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

89

6. Solution algorithm In this section, we discuss solution approaches that we could use to solve the ITND model. The ITND problem is without a doubt NP-hard; it is, in fact, harder than integer programming since it involves fixed-point equation (12) as a constraint. Based on our results in Section 5, we attack the hardness from two directions: (i) In Section 6.1, we drop the integrality of the network design vector, x, in the ITND-2 model and then solve the relaxed model using a nonlinear optimization algorithm/solver—this is the same as solving a continuous ITND problem. This method, however, possibly leads to local optima due to the non-convexity of the model. (ii) In Section 6.2, we propose an exact branch-and-bound based global optimization algorithm to solve the ITND-1R model, which returns a global minimizer. The first approach can be thought of as a heuristic (since it might find a local optimal solution to the ITND-2 model, from which we can obtain a feasible but possibly suboptimal solution to the ITND model), while the second one generates a globally optimal solution to the ITND-1R model, thus, providing a global optimal solution to the ITND model according to Theorem 2. However, when developing optimization algorithms, we often have to choose with the classical dilemma: we either sacrifice solution quality for time efficiency or the other way around. The first method may return a local optimizer but requires much shorter computational time than the second one does to find a solution. We explore both the methods and would test the trade-off between solution quality and time efficiency for them. 6.1. Nonlinear optimization method We obtain the following ITND-2Rx model by dropping the integrality of x from the ITND-2 model:

(ITND-2Rx )

C2∗x := min C2 (z ) z

subject to: (10 ) − (12 ), (16 ), (23 ), x ∈ conv(X ). Remark 9. It follows from Proposition 5 that (11) is redundant once (12) is satisfied. We hereby include (11) so that we can legitimately drop the “∧” sign in (1) and (2) when solving the ITND-2Rx model. Note that the ITND-2Rx model is a non-convex constrained nonlinear optimization model. We are interested in solving it using an iterative-search type algorithm such as the trust-region algorithm, sequential quadratic programming, or interiorpoint method (Nocedal and Wright, 2006), which often requires the objective function as well as its gradient as inputs. We will take the advantage of state-of-the-art optimization techniques and solve the INTD-2Rx model utilizing a nonlinear optimization solver that implements one of the above-referenced iterative-search methods, e.g., Matlab, into which we will feed the objective function (22), gradient (25), Jacobian (26), and an initial search point. 6.1.1. Objective function and its gradient In this section, we present the compact form of the objective function, C2 (z), defined in (22) and derive first-order derivatives (gradient) for the objective function. Define ba := (bia , i ∈ Ia ) and ea := (eia , i ∈ Ia ) for each a ∈ Ac ∪ At . It follows that

C2 (z ) =



a∈Ac ∪At

(va ba ya + ea ya ) +



(xa − 1 )e1a + ρ

a∈A2



va ,

a∈Ad

from which the first-order derivatives are given by,

  ∂ C2 (z ) ba ya if a ∈ Ac ∪ At , := ρ if a ∈ Ad , ∂va  ∂ C2 (z ) e1a if a ∈ A2 , := 0 if a ∈ A0 ∪ A1 , ∂ xa ∂ C2 (z ) := va ba + ea , ∀ a ∈ Ac ∪ At . ∂ ya

(25)

6.1.2. Jacobian of the fixed-point equation We derive in this section the Jacobian of equality constraint (12). Let Jx f(x ) := [∂ fi (x )/∂ x j ] ∈ Rm×n denote the Jacobian of any vector-valued function f := ( f1 , f2 , . . . , fm ) ∈ Rm with respect to (w.r.t.) x ∈ Rn , where the (i, j)-entry, ∂ fi (x)/∂ xj , is the derivative of fi (x) w.r.t. xj .

90

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Under Assumption 2, qw is differentiable for all w ∈ W. For each w ∈ W, let Jg qw (gw ) denote the Jacobian of qw w.r.t. gw , let Jv u(v, x) denote the Jacobian of u w.r.t. v, and let Jx u(v, x) denote the Jacobian of u w.r.t. x. It follows that the Jacobian matrix of the fixed-point operator F is obtained as

Jv,x F (v, x ) = [Jv F (v, x ) − I, Jx F (v, x )],

(26)

where I ∈ R|A|×|A| is the identity matrix,

Jv F (v, x ) :=



Dw w Jg qw (gw (v, x )) w Jv u (v, x ),

w∈W

Jx F (v, x ) :=



Dw w Jg qw (gw (v, x )) w Jx u (v, x ).

w∈W

The derivation of Jg qw (gw ) is provided in Appendix A for the case where the choice operator qw defined in (5) is described by a nested logit model (Ben-Akiva and Lerman, 1985; Train, 2003). 6.1.3. Determining an initial point In this section,  we discuss  how to determine a fairly good initial solution to the ITND-Rx model. Define y0 := yi0,a , i ∈ Ia . Let z0 = (v0 , x0 , y0 , vi0,a , i ∈ Ia , a ∈ Ac ∪ At ) denote the initial point fed into the solver, where x0 is the initial network design vector and v0 is the initial network flow distribution. It follows from Theorem 4 that the ITND-2R model provides a lower bound on the optimal objective value of the ITND-2 model. The network design vector, as part of the optimal solution to the ITND-2R model, is a good candidate that we could use as the initial network design vector for solving the ITND-2Rx model. Thus, we first solve the ITND-2R model to obtain x0 from the optimal solution and then find v0 by solving the following system of nonlinear equations:

F ( v0 , x 0 ) − v0

=

0.

(27)

It follows from Proposition 3 that v0 is the unique solution that satisfies the above equation. Again, we solve (27) using a solver (e.g., Matlab). Let yi0,a = 1 and vi0,a = v0,a if v0,a ∈ [ai−1 , ai ); let yi0,a = vi0,a = 0 otherwise for i ∈ Ia . We have thus far found an initial feasible point z0 for solving the ITND-2Rx model. We present the nonlinear optimization algorithm in Algorithm 1 . Algorithm 1: Nonlinear optimization algorithm.





Input: Ia , bia , eia , ai , i ∈ {1, 2, . . . , Ia }, Ba , a , γa , pa , ta0 , θa0 for all a ∈ Ac ∪ At , gw,0w and Dw for w ∈ W, (Bh , h ∈ H ), α , β , B , D , ρ , , F ( v , x ) . Output: Optimal network design vector x. Step 0. Solving the ITND-2R model. Solve the ITND-2R model and obtain the optimal network design vector x0 using an MILP solver. Step 1. Determining an initial point. Find an initial solution, z0 , by solving the fixed-point equation (27) with x0 obtained in Step 0. Step 2. Solve the ITND-2Rx model. Solve the ITND-2Rx model using a solver with the inputs of objective function (22), gradient (25), and Jacobian (26) of equality constraint (12) and find an optimal solution z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ), in which x is the desired optimal network design vector.

6.2. A branch-and-bound based global optimization algorithm In this section, we present a branch-and-bound based global optimization algorithm to solve the INTD-1R model, by which we can obtain an (integer) optimal solution to the ITND model according to Theorem 2. Denote by  := |A1 | + |A2 | the total number of decision variables in (xa , a ∈ A1 ∪ A2 ). We order all the entries of vector (xa , a ∈ A1 ∪ A2 ) in a prescribed sequence and index all its entries as (x1 , x2 , . . . , x ), where xi ∈ {0, 1} for all i ∈ {1, 2, . . . , }. Thus, the entries of the network design vector are indexed as x := (x1 , x2 , . . . , x , xa = 1, a ∈ A0 ). We consider a binary tree generated based on values of vector (x1 , x2 , . . . , x ). Each node on the tree represents a specified/determined value, 0 or 1, for an entry in (x1 , x2 , . . . , x ). Let a dummy node be the tree’s unique root. Each path from the root to a bottom leaf of the tree corresponds to a complete value for vector (x1 , x2 , . . . , x ). Thus, the tree’s hight is  and there are (2+1 − 1 ) nodes on the tree (including the root). Denoted by ν ∈ {0, 1, . . . , 2+1 − 2} a node on the tree, where ν = 0 denotes the root. Each node has a depth ν ∈ {0, 1, 2, . . . , }. We will implement a depth-first search strategy when exploring the nodes on the binary tree. After a node is explored/visited, the node will eventually be branched into two children nodes if necessary and applicable, one representing 0 and the other representing 1. We will then examine its two children nodes. Let B ⊂ {0, 1, . . . , 2+1 − 2} denote the set of

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

91

nodes that have been visited, can be further branched (i.e., depth is less than ), but not yet branched or pruned. Let C denote the set of leaf nodes that cannot further be branched and represents a complete value for (x1 , x2 , . . . , x ). For each node ν ∈ B with depth v , we have a path from the root to the node. Along the path, each entry from x1 to xν has been specified a value of 0 or 1. Let χν ∈ {0, 1}ν denote the vector of specified values along the path, which fixes entry values for part of the network design vector x. For each ν ∈ B, we define Sν as the set of complete values of network design vector x, in which the first ν values have already been specified, namely, Sν := {x ∈ X : (x1 , x2 , . . . , xν ) = χν }. We present the branch-and-bound based global optimization algorithm in Algorithm 2 . Algorithm 2: Branch-and-bound based global optimization algorithm.





Input: Ia , bia , eia , ai , i ∈ {1, 2, . . . , Ia }, Ba , a , γa , pa , ta0 , θa0 for all a ∈ Ac ∪ At , gw,0w and Dw for w ∈ W, (Bh , h ∈ H ), α , β , B, D, , F (v, x ), ρ , , initial set C = B = {0}, ub = ∞, lb = 0, tolerance error ε > 0. Output: Optimal solution z∗ to the ITND-1R model. Step 0. Determining the initial point. Solve the ITND-2R model to obtain x0 and then findv0 that satisfies (27) using a solver. Create z0 as is described in Section 6.1.3. Let ν0 denote the node associated with x0 . Set upper bound ubν0 ← C˜1 (z0 ), C ← C ∪ {ν0 }, and ub ← min{ubν0 , ub}. Step 1. Selecting a node. Choose ν ∈ arg maxν  ∈B ν  such that ν has the largest depth among all the nodes in B. Step 2. Branching. Breed two children nodes from the parent node ν . The left children node is denoted by ν1 with value 0 and the right one is represented by ν2 with value 1. Update χν1 ← (χν , 0 ) and χν2 ← (χν , 1 ). Set ν1 = ν2 ← ν + 1. Set B ← B \ {ν}. Step 3. Validating constraints. Execute the following loop and then go to Step 6.for j = 1 to 2 do ν˜ ← ν j ; Form a complete solution xν˜ ∈ Sν˜ , where x = 0 for all  ∈ {ν˜ + 1, ν˜ + 2, . . . , }; if xν˜ satisfies (10) then If ν˜ < , go to Step 4; go to Step 5 otherwise; else Prune ν˜ ; Step 4. Calculating the lower bound. Solve the ITND-2R model with (21) replaced by x ∈ Sν˜ using a mixed-integer linear optimization solver and set lbν˜ as the optimal objective value. If lbν˜ < ub, then B ← B ∪ {ν} ˜ ;prune ν˜ otherwise. Step 5. Updating the upper bound. Find v such that F (v, xν˜ ) = v using a solver.Let yia = 1 and via = va if va ∈ [ai−1 , ai ); let yia = via = 0 otherwise for i ∈ {1, 2, . . . , Ia }. Set zν˜ ← (v, xν˜ , yia , via , i ∈ Ia , a ∈ Ac ∪ At ). Set upper bound ubν˜ ← C˜1 (zν˜ ), C ← C ∪ {ν} ˜ , and ub ← min{ub, ubν˜ }. Step 6. Stopping test. Set lb ← minν  ∈B lbν  . If B = ∅ or |ub − lb| < ε , find z∗ ← zν ∗ where ν ∗ ∈ arg minν  ∈C ubν  , and stop; go to Step 1 otherwise.

6.3. Heuristic network design policies From solving the ITND-2Rx model using Algorithm 1, we could possibly obtain a fractional optimal network design vector x. If the network planner is keen to solve a continuous ITND problem, we have in fact accomplished the task—xa is explained as the optimal fraction of γ a that the network planner should add to the existing capacity of link a. One may wonder how we could use an optimal but fractional network design vector x if the network planner insists on solving a discrete ITND problem. In this section, we discuss two simple heuristic policies—ranking-based and sampling-based policies—which transform a fractional x into an integral network design vector which the network planner can put into practical use. If any of the entries in x is integral, then use it. We focus on transforming fractional entries into binary values with available budget, which is the total budget B less the sum of the costs for establishing/improving all a with xa = 1. Thus, we assume, without loss of generality, that all entries of x are fractional in this section. Ranking-based heuristic policy. Let xπ denote the integral network design vector transformed from x = (xa , a ∈ A0 , xa , a ∈ A1 ∪ A2 ) using the ranking-based heuristic policy described in Algorithm 3 . Sampling-based heuristic policy. We present a sampling-based heuristic policy in Algorithm 4 . Let xϖ denote the integral network design vector transformed from a fractional x. 7. Numerical experiments In this section, we offer two numerical experiments to test our solution algorithms, in which one experiment is carried out with a small-size intermodal network, while the other is provided for an intermodal network of a practical size. For both the networks, we will implement Algorithms 1 and 2. After we obtain the optimal (possibly fractional) network design vector x using Algorithm 1, we will transform it into an integral design vector xπ using Algorithm 3 and an integral

92

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Algorithm 3: Ranking-based heuristic policy. Input: (Ba , a ∈ A1 ∪ A2 ), x obtained from Algorithm 1, and B. Output: xπ . Step 0. Setting variables for A0 . Set Xa ← 1 for all a ∈ A0 . Step 1. Ranking network design variables. Re-arrange the entries of (xa , a ∈ A1 ∪ A2 ) according to a descending order. Let (xa(1 ) , xa(2 ) , . . . , xa( ) ) denote the vector of ordered network design variables, where xa(1 ) ≥ xa(2 ) ≥ . . . ≥ xa( ) , a( j ) ∈ A1 ∪ A2 is the index of link a ∈ A1 ∪ A2 with xa being the j-th largest amongst all the entries in (xa , a ∈ A1 ∪ A2 ), and j ∈ {1, 2, . . . , }.  Step 2. Checking up the budget constraint. If a∈A1 ∪A2 Ba ≤ B, set Xa ← 1 for all a ∈ A1 ∪ A2 , set xπ ← (Xa , a ∈ A ), and   j stop; otherwise find ∗ ← sup j ∈ {1, 2, . . . , } : j =1 Ba( j ) ≤ B . Step 3. Setting integer values. Set Xa( j ) ← 1 for all j ∈ {1, 2, . . . , ∗ } and Xa( j ) ← 0 for all j ∈ {∗ + 1, . . . , }. Set xπ ← (Xa , a ∈ A ) and stop.

Algorithm 4: Sampling-based heuristic policy. Input: (Ba , a ∈ A1 ∪ A2 ), x obtained from Algorithm 1, and B. Output: x . Step 0. Setting variables for A0 . Set Xa ← 1 for all a ∈ A0 . Step 1. Sampling process. For each a ∈ A1 ∪ A2 , generate a Bernoulli data point Xa ∼ Bernoulli (xa ) such that P(Xa = 1 ) = xa . Set A+ ← {a ∈ A1 ∪ A2 : Xa = 1} and + ← |A+ |.  Step 2. Checking up the budget constraint. If a∈A+ Ba ≤ B, set x ← (Xa , a ∈ A ) and stop; otherwise go to Step 3. Step 3. Ranking the network design variables. Re-arrange the entries of (xa , a ∈ A+ ) according to a descending order. Let (xa(1 ) , xa(2 ) , . . . , xa(+ ) ) denote the vector of ordered network design variables, where xa(1 ) ≥ xa(2 ) ≥ . . . ≥ xa(+ ) , a( j ) ∈ A+ is the index of link a ∈ A+ with xa being the j-th largest amongst all the entries of A+ , and j ∈ {1, 2, . . . , + }.   j Step 4. Determining positive variables. Find ∗+ ← sup j ∈ {1, 2, . . . , + } : j =1 Ba( j ) ≤ B and set Xa( j ) ← 1 for all j ∈ {1, 2, . . . , ∗+ } and Xa( j ) ← 0 for all j ∈ {∗+ + 1, . . . , + }. Step 5. Setting integer values. Set x ← (Xa , a ∈ A ) and stop.

design vector xϖ using Algorithm 4. Let x∗ denote the globally optimal network design vector that we obtain from solving the INTD-1R model using the branch-and-bound based global optimization Algorithm 2. Then, we need to compare the solution qualities for a set of three solutions:  := {xπ , xϖ , x∗ }. We suggest the following work flow to compare the three solutions: (i) For x∗ , we first solve F (v∗ , x∗ ) − v∗ = 0 to obtain v∗ and then compute the objective value C∗ := C(v∗ , x∗ ) in the ITND model. (ii) For xπ , we first solve F (vπ , xπ ) − vπ = 0 to obtain vπ and then compute the objective value Cπ := C(vπ , xπ ). (iii) To dampen the effect of random samples on the quality of xϖ , we compute the sample average approximation (SAA) of the objective value C(v, x). First, we generate a sample of N independent and identically distributed (i.i.d.) data points, {Xaj }Nj=1 for each a ∈ A1 ∪ A2 , where each Xaj ∼ Bernoulli(xa ) for all j ∈ {1, 2, . . . , N}. Second, for each j ∈ {1, 2, . . . , N}, implement Algorithm 4 to obtain an integral network design vector x . Third, solve the fixed-point equation j F ( v , x ) − x = 0 to obtain v for each j. Lastly, we compute the SAA j j j j

C  :=

N 1  C ( v j , x j ), N j=1

and will consider it as the total cost under xϖ . We set the sample size N = 10, 0 0 0 in our experiments. Both Algorithms 1 and 2 are coded with Matlab 2015a and executed under 64-bit Windows 7 on a PC with Intel i5 2.60 GHz CPU and 8 GB RAM. We choose the sequential quadratic programming algorithm embedded in fmincon of Matlab to solve the ITND-2Rx model, use Gurobi 6.5.0 with default settings as the MILP solver to solve the ITND-2R model, and utilize the fsolve function in Matlab to solve the fixed-point equation (27). 7.1. Nested logit model for intermodal operators’ route choice We consider that the choice operator q = (qrw , w ∈ W, r ∈ Rw ) defined in (5) is given by a nested logit model (Train, 2003). We provide theoretical foundations and assumptions for defining the nested logit model in Appendix A. In our nu-

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

93

Fig. 3. The intermodal network for numerical tests under Case I. Table 1 Parameter values for link capacities, construction budgets, and utilities under Case I.

Spoke carrier link Hub carrier link Existing transfer link Potential transfer link

Ba

pa

ta0

θa0

a

γa

0 0 0 1

0.5 0.3 1 1

0.1 0.2 – –

1 1 0.5 0.5

15 20 10 0

0 0 0 10

merical tests, we consider that the set of routes between each O–D pair w ∈ W, i.e., Rw , is partitioned into k = 2 nests, R1w = {0w } and R2w = Rw \ {0w }. 7.2. Case I In this case, we consider a small-size intermodal network shown in Fig. 3, in which the number along each link shows the link’s id and the number inside a circle/rectangle indicates the node’s id. It follows that the set of O–D pairs, W = {(1, 5 ), (1, 6 ), (2, 5 ), (2, 6 )}, the set of carrier links, Ac = {1, 2, . . . , 10}, and the set of no-choice links Ad = {11, 12, 13, 14}. We have that At = A(3 ) ∪ A(4 ), where, A(3 ) = {(1, 5 ), (1, 7 ), (1, 8 ), (3, 5 ), (3, 7 ), (3, 8 ), (6, 5 ), (6, 7 ), (6, 8 )} and A(4 ) = {(2, 6 ), (2, 9 ), (2, 10 ), (4, 6 ), (4, 9 ), (4, 10 ), (5, 6 ), (5, 9 ), (5, 10 )}. Note that each number inside the parenthesis represents a carrier link’s id. We consider that A0 = Ac ∪ (A(3 ) \ {(6, 8 )} ), A1 = ∅, and A2 = A(4 ) ∪ {(6, 8 )}. Thus,  = 10. 7.2.1. Data  We set demands D(1,5 ) = 15, D(1,6 ) = 3, D(2,5 ) = 5, D(2,6 ) = 7. Let Ia = 3, ba = (2, 1, 0.5 ), and ea = 0,

10 20 3 , 3



for all a ∈

Ac ∪ At . We use the total budget B = 6, α = 1, β = 1.2, and ρ = 100. Set τw = 0.1, τw1 = τw2 = 0.2, and gw,0w = −150.36 for all w ∈ W. Set tolerance error ε = 10−7 . Table 1 shows the parameter values for link capacities, construction budgets, and utilities under Case I. 7.2.2. Computational results Table 2 shows partial computational results for transfer links that are to be established. It takes less than 4 s for Algorithm 1 to terminate, while we spend 26.4 s before Algorithm 2 stops and returns a global optimal solution. Fig. 4 shows the comparison of C∗ , Cπ , and Cϖ under Case I. 7.3. Case II For this case, we implement our numerical tests on an intermodal network with three models of carrier links: rail, truck, and sea, as shown in Fig. 5. The network consists of 35 nodes (with node ids next to them), out of which there are three types of intermodal hubs: (i) inland intermodal hubs such as rail-rail or rail-truck dry ports, (ii) seaports, and (iii) border hubs between geographically neighbored countries. Table C.5 lists the id, head, tail, mode, and length of each of the 48 carrier links in the network. Table C.6 shows the demands for all O–D pairs under consideration.

94

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104 Table 2 Computational results for transfer links under Case I. Algorithm 1

Algorithm 2

a ∈ A2

xa

va

y1a

(6,8) (2,6) (2,9) (2,10) (4,6) (4,9) (4,10) (5,6) (5,9) (5,10)

0.48 0.00 1.00 0.68 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 7.68 1.50 0.00 0.00 0.00 0.00 0.00 0.00

1 1 0 1 1 1 1 1 1 1

y2a

y3a

xa

va

0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0

1 0 1 0 0 0 0 0 0 0

0.00 0.00 7.74 0.00 0.00 0.00 0.00 0.00 0.00 0.00

4s

Time

26.4 s

Fig. 4. Comparison of C∗ , Cπ , and Cϖ under Case I. Table 3 Parameter values for link capacities, construction budgets, and utilities under Case II.

Carrier link Carrier link Carrier link Existing transfer link Potential transfer link

Mode

Ba

s0a

p¯ a

θa0

a

γa

Truck Rail Sea – –

0 0 0 0 1

27.7 33.3 21.6 – –

0.8 0.1 0.2 27.0 27.0

1 1 1 12 12

800 800 800 300 0

0 0 0 0 300

Let H˜ := {2, 6, 8, 9, 11, 13, 15, 18, 22}. We consider that A2 = ∪h∈H˜ A(h ) and that all other links belong to A0 . For simplicity, we ignore cargo transfer processes at all other intermodal hubs not in H˜ . We have that  = 19 transfer links to be established at all hubs in H˜ . 7.3.1. Data Let Ia = 3, ba = (2, 1, 0.5 ), and ea = (0, 40, 80 ) for all a ∈ Ac ∪ At . We use the total budget B = 15, α = 1, β = 1.2, and ρ = 10 0 0. Set τw = 0.1, τw1 = τw2 = 0.2, and gw,0w = −3 × 105 for all w ∈ W. Set the tolerance error to be ε = 10−7 . Table 3 shows the parameter values for link capacities, construction budgets, and utilities, in which s0a is the average transportation speed on link a ∈ Ac , p¯ a denotes the transportation cost per unit distance over link a if a ∈ Ac , and it represents the transfer cost if a ∈ At . 7.3.2. Computational results Table 4 shows partial computational results obtained using Algorithms 1 and 2 for all transfer links to be established under Case II. Fig. 6 shows the comparison of C∗ , Cπ , and Cϖ under Case II. Under Case II, it takes 295 s for Algorithm 1 to stop, while Algorithm 2 runs for 1434 s before termination.

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

95

Fig. 5. The intermodal network for numerical tests under Case II.

7.4. Findings and discussions We discuss in this section what we have learned from our numerical experiments. We obtain the following interesting findings: (i) The first observation that we made from Tables 2 and 4 is that only one of the three segment variables, i.e., y1a , y2a , and y3a , is one and the other two are zeros, which is consistent with Proposition 6. In addition, for all those links with xa = 0, we have that the cargo flows loaded on them are also zero, which is consistent with Proposition 4. (ii) Second, we learn from Figs. 4 and 6 that C∗ < Cπ and C∗ < Cϖ , which agrees with the global optimality of x∗ . (iii) Third, a clear trade-off between solution quality and computational time shows up in our numerical tests, as is often observed in the numerical studies on optimization algorithms. Under Case II, it seems to us that the network design vector xπ (obtained based on the ranking-based heuristic policy π ∗ and Algorithm 1) provides quite a good solution that leads to a total cost C C−C = 0.8% higher (worse) than the glob∗ −295 ally minimum cost. However, being only “0.8% worse”, Algorithm 1 saves us 1434 = 3.86 times the computational 295 time compared with the global optimization Algorithm 2.

96

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104 Table 4 Computational results for transfer links under Case II. Algorithm 1

Algorithm 2

a ∈ A2

xa

va

y1a

(1,2) (1,3) (4,8) (10,21) (11,21) (14,21) (21,23) (22,23) (12,15) (13,15) (8,27) (23,27) (27,28) (26,36) (28,36) (29,36) (32,36) (34,35) (41,35)

1.00 0.00 0.00 0.00 0.10 0.03 0.13 0.00 1.00 1.00 0.00 0.13 0.13 0.38 0.13 1.00 0.09 0.29 0.05

190.38 0.00 0.00 0.00 30.83 7.88 38.71 0.00 57.38 12.31 0.00 38.71 38.71 114.06 38.71 300.00 26.94 86.21 15.74

0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 0 1 1 1

Time

y2a

y3a

xa

va

0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0

1 0 0 0 1 1 1 0 1 1 0 1 1 0 1 1 0 1 0

190.38 0.00 0.00 0.00 31.53 8.07 39.60 0.00 56.67 12.12 0.00 39.60 39.60 0.00 39.60 300.00 0.00 168.33 0.00

295 s

1434 s

Fig. 6. Comparison of C∗ , Cπ , and Cϖ under Case II. .

Under Case I, the network design vector xϖ (obtained from the sampling-based heuristic policy and Algorithm 1) gives  ∗ a total cost that is C C−C = 1.04% worse (higher) than the optimal cost. Again, it requires a much shorter time to run ∗ Algorithm 1 (4 s) than the time to execute Algorithm 2 (26.4 s), as expected. At this time point, we feel compelled to arrive at a conclusion that, though Algorithm 1 might return a locally optimal solution, we could obtain quite a good heuristic solution to the INTD model using either the sampling-based or rankingbased heuristic policy, depending on the size of the problem. On one hand, we might only lose insignificant amount of money using the nonlinear optimization Algorithm 1; on the other hand, what excites us most is that we would save a considerable amount of computational time to reach a reasonably good solution. Note that the INTD problem is a planning-level problem and it can be solved offline with plenty of hours using powerful computational resources. Both Algorithm 1 and Algorithm 2 can find good applications when solving the INTD problem. With all the results and discussions above, we might be in a good position to offer suggestions on how to choose a proper algorithm to solve the ITND model: the global optimization Algorithm 2 should be used if one would like to sacrifice time for a better solution quality and the nonlinear optimization Algorithm 1 seems appropriate if one is inclined to embrace a fast, though slightly inferior, solution.

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

97

8. Concluding comments In this paper, we studied the discrete intermodal freight transportation problem with three prominent features: economies of scale, congestion effects, and route choice behaviors of intermodal operators. The problem was first formulated as a mixed-integer nonlinear and non-convex programming model involving both discrete and continuous variables, piecewise linear cost functions, and a fixed-point constraint, making it a hard problem. To solve the model, we formulated a number of equivalent and relaxed models with stepwise removal of constraints and the integrality of integer decision variables. The equivalences between models were further justified by theoretical results. This does not only ease our problem but also provides insights into how to simplify a minimization model with piecewise linear functions as a component of its objective function. We offered two solution strategies: a nonlinear optimization algorithm with policies generating heuristic solutions to the ITND problem and a global optimization algorithm. Several managerial insights are gained. In particular, we found that it would not lose much for the network planner to use the nonlinear optimization algorithm together with the sampling-based heuristic policy to solve the ITND problem of a mall size. We also learned that it gives a quite reasonable solution within a short time by solving the ITND problem of a larger size using the nonlinear optimization algorithm with the rankingbased heuristic policy. Furthermore, if the network planner has plenty of time and is a quality lover, it would be the best to run the global optimization algorithm offline. For future studies, we would be interested in applying our methodology to maximize the market share of a logistics hub by designing its hinterland network and discussing approximation methods to solve the non-convex ITND model. Acknowledgments We are indebted to the Editor-in-Chief and our two anonymous referees for their thoughtful comments that have helped substantially improve this work. The authors thank Shuaian Wang and Hua Wang for their valuable input on this work. This study is supported by the research project “Liner Shipping Container Slot Booking Patterns and their Applications to the Shipping Revenue Management” funded by the NOL Fellowship Programme, Singapore. Appendix A. Nested logit model We provide details of the nested logit model in this section. Interested readers are referred to Train (2003) for more theories and discussions of discrete choice models, including the nested logit model. For each O–D pair w ∈ W, we partition Rw into Kw + 1 nests/groups, R0w , R1w , R2w , . . . , RKww , where Kw ∈ N, R0w = {0w } and

Rw =

Kw

Rkw and Riw ∩ Rwj = ∅ for i = j,

k=0

where the nests are classified according to similar but unobservable attributes of the routes. The partition can be obtained according to some principles, e.g., one principle could be that we distinguish the routes that do not contain ocean legs and those traversing through ocean legs. We could also classify an intermodal route into a nest according to whether or not the k = 1 if r ∈ Rk ; zero otherwise. route traverses a port. For w ∈ W, k = 0, 1, . . . , Kw and r ∈ Rw , let δw,r w k For each w ∈ W, let ξw be the random variable that represents the nest-related random error in the systematic utilities k of routes in Rkw and ξw,r be the random variable that denotes the alternative-related random error for r ∈ Rkw . We need Assumption 3 to formulate the nested logit model. Assumption 3. The following is assumed. k are independent for all w ∈ W, k ∈ {0, 1, . . . , K } and r ∈ Rk and independent of g . (1) ξwk and ξw,r w wr w k , r ∈ R } is a sequence of i.i.d. Gumbel random variables supported on k ⊂ R with scale parameter τ k ∈ (0, ∞ ) (2) {ξw,r w w,r w and location 0 for all w ∈ W and k = 0, 1, . . . , Kw . K K w w (3) {ξwk }k=1 is distributed such that {maxr∈Rk g˜w,r }k=0 is a sequence of i.i.d. Gumbel random variables with scale paramew

ter 0 < τw < min{τwk : k = 0, 1, . . . , Kw } for all w ∈ W. Each ξwk is supported on kw ⊂ R for w ∈ W and k ∈ {0, 1, . . . , Kw }.

k , w ∈ W, k ∈ {0, 1, . . . , K }, r ∈ Rk ). For any w ∈ W, the random route utility g ˜w,r is given by Define ξ := (ξwk , ξwr w w

g˜w,r (v, x, ξ )

:=

gw,r (v, x ) +

Kw 

k k δw,r (ξw,r + ξwk ) = gw,r (v, x ) + ξw,r , ∀w ∈ W, r ∈ Rw .

k=1

By the nested logit model, the probability of intermodal operators choosing r ∈ Rw , where w ∈ W, defined in (5), is given by

qrw (gw (v, x ))

=

Kw  k=0

k δw,r qrw,k (gw (v, x ))qw,k (gw (v, x )),

98

where

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

 τwk gw,r (v, x )  k , (gw (v, x )) :=  i∈Rkw exp τw gw,i (v, x )   exp τw g¯ w,k (v, x )  , qw,k (gw (v, x )) := K w ¯ j=0 exp τw gw, j (v, x ) exp

qrw,k



subject to the understanding that 0/0 = 0, where



g¯ w,k (v, x ) :=

1

τwk



ln

exp



τ

k w gw,r

( v, x )

 

.

r∈Rkw

Proposition 7. The Jacobian of qw w.r.t. gw , Jg qw (gw ), for each w ∈ W, is derived as,

⎧ −τw qrw qsw ⎨ ∂ qrw r s q q (τw (1 − qw,k ) − τwk ) = ∂ gw,s ⎩ r wr w,k τw qw qw,k (1 − qwk ) + τwk qrw (1 − qrw,k )

if if if

r ∈ Rkw , s ∈ Rlw , k, l ∈ {0, 1, . . . , Kw }, k = l, r, s ∈ Rkw , k ∈ {0, 1, . . . , Kw }, r = s, r = s ∈ Rkw , k ∈ {0, 1, . . . , Kw }.

Proof. We suppress w for notational simplicity. If r ∈ Rk , s ∈ Rl , k, l ∈ {0, 1, . . . , K }, k = l, it follows that

  exp τ k gr ∂ qr exp (τ g¯ k )   K   =  kg ∂ gs exp τ k i i∈R j=0 exp τ g¯ j    l  l exp τ k gr −τ exp (τ g¯ k + τ g¯ l ) 1 τ exp τ gs      =   2 τ l  k l K i∈Rk exp τ gi i∈Rl exp τ gi j=0 exp τ g¯ j = −τ qr qs .

If r, s ∈ Rk and r = s, it follows that

  −τ k exp τ k gr + τ k gs ∂ qr exp (τ g¯ k )   =   2 K ∂ gs k j=0 exp τ g¯ j k exp τ g i∈R

i





exp(τ k gs )

τ exp(τ g¯ k )  exp τ k g τ k gr ( i) i∈Rk   + k i∈Rk exp τ gi exp

 ( τ k gs ) τ g¯ j − τ exp(τ g¯ k + τ g¯ k )  expexp ( τ k gi ) i∈Rk  K  2 j=0 exp τ g¯ j

K

j=0

exp



= qr qsk (τ (1 − qk ) − τ k ). If r = s ∈ Rk , it follows that

  exp τ k gr ∂ qr exp (τ g¯ k )   K   =  k ∂ gr i∈Rk exp τ gi j=0 exp τ g¯ j    k   k   k  exp τ k gr i∈Rk exp τ gi − exp τ gr exp τ gr k = τ qk   2 k i∈Rk exp τ gi    exp(τ g¯ k )qrk Kj=0 exp τ g¯ j − exp(τ g¯ k ) exp(τ g¯ k )qrk r + τ qk  K  2 j=0 exp τ g¯ j τ k qk (qrk − (qrk )2 ) + τ qrk (qr − (qk )2 qrk ) = τ qr qrk (1 − qk ) + τ k qr (1 − qrk ), =

which proves the result.



Appendix B. Technical proofs Proof of Proposition 1. Let v ∈ R|A| be such that v0w = Dw for all w ∈ W and va = 0 for all a ∈ Ac ∪ At . Thus, v ∈ v and v = ∅. Since the no-choice link 0w is alsoa route connecting O–D pair w ∈ W, let f0w = Dw for all w and we have f ∈  f , which implies  f = ∅. Note that |va | ≤ w∈W Dw and | fw,r | ≤ w∈W Dw . In addition, the closeness of v and f follows from the linear transformations in their definitions. Thus, both the sets are compact.  Proof of Proposition 2. By Remark 3, v × conv(X ) = ∅. Consider any (v, x ) ∈ R|A| × conv(X ). Define fw := Dw qw (gw (v, x)). Since qw (gw (v, x ))1 = 1 by Remark 4, it follows fw 1 = Dw . Define f := (fw , w ∈ W ). We have that F (v, x ) =

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104



|R |

w∈W w fw = f. Thus, there exists f ∈ R+ and completes the proof. 

99

such that F (v, x ) = f, fw 1 = 1 for all w ∈ W, which implies F(v, x) ∈ v

¯ |Rw | , we have from Cascetta (2009, pp. 135–137) that, if g1 1 , Proof of Lemma 1. Consider any w ∈ W. For any g1 , g2 ∈ R 2 1 2  1 g 2 < ∞ for all r ∈ Rw , then (qw (g ) − qw (g )) (g − g2 ) ≥ 0. Note that it follows from Remark 5 and our notational conventions that the above inequality holds if gir = −∞ for some r ∈ Rw and i ∈ {1, 2}. This completes the proof.  Proof of Proposition 3. It follows from Remark 2 that conv(X ) = ∅. Choose any x ∈ conv(X ). We first show the existence.  Trivially, v is a nonempty convex and closed set by Proposition 1. For each v ∈ v , it follows v2 ≤ w∈W Dw . Thus, v is a nonempty convex and compact set. The continuity of F( ·, x) follows from the continuities of ta (·, xa ), a ∈ A and the choice function qrw (· ) for all r ∈ Rw according to Assumption 2. By the Brouwer fixed point theorem, there exists a v ∈ v such that v = F (v, x ). We next show the uniqueness. Assume by contradiction that, given x ∈ conv(X ), there exist v1 , v2 ∈ v such that v1 = v2 and vi = F (vi , x ) for i = 1, 2. For any w ∈ W, by Lemma 1, it obtains

(qw (gw (v1 , x )) − qw (gw (v2 , x ))) (gw (v1 , x ) − gw (v2 , x )) ≥ 0, which yields

(Dw qw (gw (v1 , x )) − Dw qw (gw (v2 , x )) (w u(v2 , x ))) − w u(v2 , x )) ≥ 0 ⇒ (Dw w qw (gw (v1 , x )) − Dw w qw (gw (v2 , x ))) [u(v1 , x ) − u(v2 , x )] ≥ 0 ⇒



(Dw w qw (gw (v1 , x )) − Dw w qw (gw (v2 , x ))) [u(v1 , x ) − u(v2 , x )] ≥ 0 ⇒

w∈W

( v 1 − v 2 )  [ u ( v 1 , x ) − u ( v 2 , x )] ≥ 0 ,

(B.1)

F ( vi , x )

vi

where the last inequality follows from = for i = 1, 2. By Assumption 1 and Eq. (3), u( ·, x) is monotone nonincreasing, which implies

( v 1 − v 2 )  [ u ( v 1 , x ) − u ( v 2 , x )] ≤ 0 , which combining (B.1) gives that

( v 1 − v 2 )  [ u ( v 1 , x ) − u ( v 2 , x )] =



(v1a − v2a )(ua (v1a , xa ) − ua (v2a , xa )) = 0.

a∈A

Since (v1a − v2a )(ua (v1a , xa ) − ua (v2a , xa )) ≤ 0 for all a ∈ A, we have that

(v1a − v2a )(ua (v1a , xa ) − ua (v2a , xa )) = 0, ∀a ∈ A. (B.2) 1 2 1 2 For a ∈ A such that va = va , we have ua (va , xa ) = ua (va , xa ) since u( ·, xa ) is defined as a function. For a ∈ A such that v1a = v2a , we also have from (B.2) that ua (v1a , xa ) = ua (v2a , xa ). Thus, u(v1 , x ) = u(v2 , x ). It follows from (7) that F ( v1 , x ) =



Dw w qw (w u(v1 , x )) =

w∈W



Dw w qw (w u(v2 , x )) = F (v2 , x ),

w∈W

which implies v1 = v2 and contradicts v1 = v2 . This proves the uniqueness and completes the proof.



Proof of Proposition 4. Consider any a ∈ A2 such that xa = 0 in the solution (v, x ) ∈ v × conv(X ). Then, if a ∈ Ac is a carrier link, then ta = ta0 /x0 = ∞. If a ∈ At is a transfer link, then a (x0 ) = a + xa γa = 0 and it follows from (2) that ta (va , xa ) = ta0 /(−va )+ = ta0 /0 = ∞. Then, ua (va , xa ) = −∞ holds. a,r For all r ∈ Rw and w ∈ W such that δw = 1, it follows from (4) that gw,r (v, x ) = −∞ and follows from Remark 5 that r qw (gw (v, x )) = 0. We have from the route-version fixed-point equation (6) that fw,r = Dw qrw (gw (v, x )) = 0, which shows (1). It follows from (7) that

va =



Dw

w∈W

which proves (2).



r∈Rw

δwa,r qrw (gw (v, x )) =

 

δwa,r fw,r = 0,

w∈W r∈Rw



Proof of Lemma 2. (1) If x = 0, we have y f (0 ) > 0 = x f (y ). Consider now x > 0. Since x < y, let λ = x/y and we have 0 < λ < 1. Thus, it follows from the concavity of f and f(0) > 0 that

f (x ) f (λy ) f (λy + (1 − λ )0 ) = = ≥ x λy λy

λ f (y ) + (1 − λ ) f (0 ) λ f (y ) f (y ) > = . λy λy y

(2) Now, consider x, y ∈ R+ . If x = y = 0, we have f (x ) + f (y ) = 2 f (0 ) > f (0 ) = f (x + y ). We consider x + y = 0 now. It follows from (1) that f (x + y )x/(x + y ) < f (x ) and f (x + y )y/(x + y ) < f (y ). Then,

f (x + y ) = which proves (2).

f ( x + y )x f ( x + y )y + < f ( x ) + f ( y ), x+y x+y

100

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

(3) If x = 0, it holds that y f (0 ) = 0 = f (y )x = 0. If f is a linear function passing through 0, then f (y )x = f (x )y. On the other hand, if f (y )x = f (x )y holds for all 0 ≤ x < y, we have two cases: (i) x = 0 (ii) x = 0 and f (y )/y = f (x )/x, which implies that f is a line segment connecting 0, x, and y. Consider now x > 0 and f not being a linear function. Since 0 < x < y, let λ = x/y and we have 0 < λ < 1. Thus, it follows from the concavity of f and f (0 ) = 0 that

f (x ) f (λy ) f (λy + (1 − λ )0 ) = = ≥ x λy λy

λ f (y ) + (1 − λ ) f (0 ) f (y ) = , λy y

where it follows from Jensen’s inequality that the inequality holds as equality if and only if y = 0, which is not true, or f is a linear function, which is not true either. Thus, f(x)y > f(y)x. (4) If xy = 0, the result trivially holds. If f is a linear function passing through 0, then f (x + y ) = f (x ) + f (y ). On the other hand, f is additive, i.e., f (x + y ) = f (x ) + f (y ) holds for all 0 ≤ x, y, we have that f (0 ) = 2 f (0 ) which implies f (0 ) = 0. We have two cases: (i) xy = 0, or (ii) it follows that

f

x + y 2

=

 

1 x+y f 2 2



+f

 x + y  2

1 f (x + y ) 2

=

=

1 ( f (x ) + f (y )), 2

which shows that f is midpoint convex. Since f is concave, it is continuous. It follows from Theorem 818 in Montesinos et al. (2015) that f is convex. Then, f is both concave and convex, and it is thus a linear function according to Exercise 1.1 in Bertsimas and Tsitsiklis (1997). Consider now xy = 0 and f not being linear. It follows from (3) that

f (x + y ) =

f ( x + y )x f ( x + y )y + < f ( x ) + f ( y ), x+y x+y

which completes the proof.



Proof of Proposition 5. Consider that (v, x) satisfies (12). Assume by contradiction that (11) is violated. Then, there exists a link a ∈ Ac ∪ At such that va > a (xa ) = a + γa xa ≥ 0, and it follows from (1) and (2) that, ta (va , xa ) = 1/(a (xa ) − va )+ = ta0 /0 = ∞. We have from (3) that ua (va , xa ) = −∞. It follows from the same argument in the proof of Proposition 4(2) that va = 0, contradicting va > 0. Thus, (11) holds.  Proof of Theorem 1. Note that Constraint (20) is redundant in the ITND-1 model. We do not consider it in our proof. First, we show that (v, x) is feasible to the ITND model. Since z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) is optimal to the ITND-1 model, then, x ∈ X , (v, x) satisfies (12), i.e., v = F (v, x ), and it follows from Proposition 2 that v = F (v, x ) ∈ v . Then, (13) is satisfied. It follows from Proposition 5 that constraint  (11) holds, which implies that (v, x) is feasible to the ITND model.   Second, we show that Ca (va ) = Iia=1 bia via + yia eia for any a ∈ Ac ∪ At . Recall that ia ∈ Ia is such that va ∈ [aia −1 , aia ). It

follows from (15)–(19) that (i) yiaa = 1, viaa = va , ya = va = 0 for all j = ia and j ∈ Ia if va = aia −1 , or (ii) yiaa −1 = 1, viaa −1 = va , j j ya = va = 0 for all j = ia − 1 if va = aia −1 . For both the cases, we have from Property 1(4) that j

Ca (va )

=

Ia  

j



bia via + yia eia .

(B.3)

i=1

Next, consider any (vˆ , xˆ ) feasible to the ITND model and we will construct a solution zˆ that is feasible to the ITND-1 model. Then, (vˆ , xˆ ) satisfies (12) and (21). For each a ∈ Ac ∪ At , there exists ıa ∈ Ia such that vˆ a ∈ [aıa −1 , aıa ). Let vˆ ıaa = vˆ a j j and yˆıaa = 1, and define vˆ a = yˆa = 0 for all j ∈ Ia and j = ıa . Denote by yˆ := (yˆia , i ∈ Ia , a ∈ Ac ∪ At ). Then, zˆ = (vˆ , xˆ , yˆ , vˆ ia , i ∈ Ia , a ∈ Ac ∪ At ) satisfies (15)–(19) and this implies that zˆ is feasible to the ITND-1 model. Due to the optimality of z and feasibility of zˆ to the ITND-1 model, it follows from (14), (B.3), and Property 1(4) that

C1 (zˆ ) =



(bıaa vˆ a + eıaa ) +

a∈Ac ∪At

=



Ca (vˆ a ) +

a∈Ac ∪At

 

a∈Ac ∪At

a∈A2

(xˆa − 1 )e1a + ρ





vˆ a

a∈Ad

vˆ a

a∈Ad

≥ C1 (z )

Ia 





bia via + yia eia +

a∈Ac ∪At i=1

=

(xˆa − 1 )e1a + ρ

a∈A2

= C (vˆ , xˆ ) =





Ca (va ) +

 a∈A2



(xa − 1 )e1a + ρ

a∈A2

(xa − 1 )e1a + ρ





va

a∈Ad

va = C ( v, x ) ,

a∈Ad

which implies C (vˆ , xˆ ) ≥ C (v, x ) and shows the optimality of (v, x) to the ITND model. Proof of Theorem 2. Choose any a ∈ Ac ∪ At . We have the following three cases:



X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Case 1.

Ia

yi = 0. Under this case, yia = via = 0 for all i. By the definition of zˆ , it holds that yˆ1a = 1, yˆia = 0 for all i = 1, i=1 a and that vˆ ia = via for all i. It follows that vˆ a = va = 0 and that

e1a + Case 2.

101

Ia  



bia via + yia (eia − e1a ) = e1a +

i=1

Ia

i=1

Ia  



bia vˆ ia + yˆia (eia − e1a ) = Ca (vˆ a ).

(B.4)

i=1

yia = 1. Under this case, let ıa be such that yıaa = 1. It follows from (15) and (20) that yia = via = 0 for all i = ıa

and va = vıaa . By the definition of zˆ , it holds that vˆ a = vˆ iaa = va , yˆiaa = 1 and yˆia = vˆ ia = 0 for all i = ia . It follows from Property 1(4) that biaa va + eiaa − e1a ≤ bıaa va + eıaa − e1a , which implies that

e1a + Case 3.

Ia

Ia  



bia via + yia (eia − e1a ) ≥ e1a +

i=1

i=1

Ia  



bia vˆ ia + yˆia (eia − e1a ) = Ca (vˆ a ).

(B.5)

i=1

yia ≥ 2. By the definition of zˆ , it holds that vˆ iaa := va , yˆiaa = 1 and yˆia = vˆ ia = 0 for all i = ia . Let Ia+ := {i ∈ Ia : k −1

yia = 1}. Let ki ∈ Ia such that via ∈ [a i

e1a +

Ia  



k , a i ). It follows from Property 1(4) of function C˜a that

bia via + yia (eia − e1a ) ≥ e1a +

i=1







bkai via + ekai − e1a = e1a +

i∈Ia+



C˜a (via ),

i∈Ia+

which with Lemma 2(4) (applied to function C˜a ) implies that,

e1a +

Ia  



bia via + yia (eia − e1a ) ≥ e1a + C˜a (va ) = e1a +

Ia  

i=1



bia vˆ ia + yˆia (eia − e1a ) = Ca (vˆ a ).

(B.6)

i=1

It follows from the three cases that zˆ = (v, x, yˆ , vˆ ia , i ∈ Ia , a ∈ Ac ∪ At ) is feasible to the ITND-1R model and it follows from (B.4)–(B.6) thatC˜1 (z ) ≥ C˜1 (zˆ ), which shows (1). Now, consider z to be optimal to the ITND-1R model. It follows from the above argument that zˆ is feasible to the ITND-1 model. Since C˜1 (z ) ≥ C˜1 (zˆ ) and z is optimal to the ITND-1R model, it follows from the feasibility of zˆ to the ITND-1R model that zˆ is optimal to the ITND-1R model. Thus, C˜1∗ = C˜1 (zˆ ). We now show that zˆ is optimal to the ITND-1 model. Consider any feasible solution, z˜ = (v˜ , x˜ , y˜ , v˜ ia , i ∈ Ia , a ∈ Ac ∪ At ), to the ITND-1 model. Thus, z˜ is feasible to the ITND-1R model. We have from the optimality of zˆ to the ITND-1R model and (B.4)–(B.6) that C˜1 (z˜ ) = C1 (z˜ ) ≥ C˜1 (zˆ ) = C1 (zˆ ). It follows from the feasibility of zˆ to the ITND-1 model that zˆ is optimal to the ITND-1 model, and therefore, C1∗ = C1 (zˆ ) = C˜1 (zˆ ). It follows from Theorem 1 that (v, x) is optimal to the ITND model. This shows (2). It follows from the above argument that C1∗ = C˜1∗ , which is the desired (3).  Proof of Proposition 6. We consider any link a ∈ Ac ∪ At . Assume by contradiction that there exists i = ia if va ∈ (aia −1 , aia ) and there exists i ∈ / {ia − 1, ia } if va = aia −1 such that 0 < yia ≤ 1. It follows from (16) that 0 ≤ yiaa < 1. Consider any δ > 0 and neighborhood Bδ (z ) := {z : z − z2 ≤ δ}. Construct z˜ := (v˜ , x˜ , y˜ , v˜ ia , i ∈ Ia , a ∈ Ac ∪ At ) such that v˜ = v, x˜ = x, y˜ia = yia − (δ /2 ) ∧ yia , y˜iaa = yiaa + (δ /2 ) ∧ yia , and y˜aj = yaj for all j ∈ {i, ia }. Since y˜ia = yia − (δ /2 ) ∧ yia ≥ yia − yia ≥ 0 and y˜iaa = yiaa + (δ /2 ) ∧ yia ≤ yiaa + yia ≤ 1, it follows that z˜ is feasible to the ITND-2 model. We have that

z˜ − z ≤



√

|y˜ia − yia |2 + |y˜iaa − yiaa |2 ≤ δ / 2,

which implies that z˜ ∈ Bδ (z ). We also have from Property 1(4), (δ /2 ) ∧ yia > 0, and va ∈ [aia −1 , aia ) that bia va + eia > biaa va + eiaa and that















C2 (z ) − C2 (z˜ ) = yia bia va + eia + yiaa biaa va + eiaa − y˜ia bia va + eia − y˜iaa biaa va + eiaa







= ((δ /2 ) ∧ yia ) bia va + eia − (biaa va + eiaa ) > 0, which contradicts the local optimality of z and proves the results.



(v, x, y˜ , v˜ ia , i

∈ Ia , a ∈ Ac ∪ At ), where y˜iaa = 1, v˜ iaa = va , and y˜ia = v˜ ia = 0 for all i Proof of Theorem 3. Consider a variable z˜ = = ia . Then, z˜ is feasible to the ITND-2 model. It follows from Corollary 2 that C2 (z ) = C2 (z˜ ), which implies that z˜ is optimal to the ITND-2 model, and therefore, C2∗ = C2 (z˜ ). Note that z˜ is feasible to the ITND-1 model. We now show the optimality of z˜ to the ITND-1 model. Consider any zˆ = (vˆ , xˆ , yˆ , vˆ ia , i ∈ Ia , a ∈ Ac ∪ At ) that is feasible to the ITND-1 model. Then, zˆ is also feasible to the ITND-2 model. Due to the optimality of z˜ to the ITND-2 model and Corollary 2 that

C2 (z˜ ) =



a∈Ac ∪At

=



a∈Ac ∪At

 



biaa va + eiaa +



biaa vˆ a + eiaa +



(xa − 1 )e1a + ρ

a∈A2



a∈A2



va = C1 (z˜ ) ≤ C2 (zˆ )

a∈Ad

(xˆa − 1 )e1a + ρ



a∈Ad

vˆ a = C1 (zˆ ),

102

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

which shows the optimality of z˜ to the ITND-1 model, and therefore, C1∗ = C1 (z˜ ) = C2 (z˜ ). It follows from Theorem 1 that (v, x) is optimal to the ITND model. It also follows from the above argument that C1∗ = C2∗ .  Proof of Theorem 4. Consider the following mathematical program:

(ITND-3 ) C3∗ := min C˜1 (z ) z

subject to: (10 ), (12 ), (15 ), (19 ), (21 ), (23 )

via ≤ yia va ∀ i ∈ {1, 2, . . . , Ia }, a ∈ Ac ∪ At .

(B.7)

The proof is carried out in four steps. Step 1. First, we show that C3∗ ≥ C2∗ . For any z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) optimal to the ITND-3 model, since it is a minimization model, it follows that (B.7) is binded, i.e., via = yia va for all link a and segment i ∈ Ia . For each link a, we construct y˜ a = (y˜ia , i ∈ Ia ) for the following two cases:  (i) If va = 0, then via = 0 for all i. Due to the optimality of z, we have Iia=1 (bia via + yia (eia − e1a ) = 0. Define y˜1a := 1 and y˜ia := 0 for all i = 1. (ii) If va > 0, define y˜ia := yia for all i, y˜ := (y˜ a , a ∈ Ac ∪ At ), and z˜ := (v, x, y˜ , via , i ∈ Ia ). Then, C˜1 (z˜ ) = C˜1 (z ) = C3∗ , which implies that z˜ is optimal to the ITND-3 model. Note that v˜ ia = y˜ia v˜ a for all a and i,  which and (15) ensure that Iia=1 y˜ia = 1. Thus, z˜ satisfies (16) and is thus feasible to the ITND-2 model. Thus, Ia 





y˜ia bia v˜ a + eia =

i=1

Ia  



bia v˜ ia + y˜ia eia = e1a +

i=1

Ia  



bia v˜ ia + y˜ia (eia − e1a ) ,

i=1

which implies that C2∗ ≤ C2 (z˜ ) = C˜1 (z˜ ) = C3∗ . Step 2. Second, we show that C3∗ ≤ C2∗ . Consider any z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) optimal to the ITND-2 model, where via := va yia for all a and i. Thus, z satisfies (15), (19), and (B.7). This implies that z is feasible to the ITND-3 model. It follows from (16) that Ia 





yia bia va + eia = e1a +

i=1

Ia  



bia via + yia (eia − e1a ) ,

i=1

which implies C2 (z ) = C˜1 (z ). Due to the feasibility of z to the ITND-3 model, we have that C2∗ = C2 (z ) = C˜1 (z ) ≥ C3∗ . We thus have that C2∗ = C3∗ . Step 3. Third, it follows from Theorems 2(3) and 3 that C2∗ = C˜1∗ . Step 4. Last, we show that the ITND-2R model is a relaxation of the ITND-3 model. For any z = (v, x, y, via , i ∈ Ia , a ∈ Ac ∪ At ) optimal to the ITND-3 model, it follows from the proof of Theorem 1 that z satisfies (11) and (24). Since z satisfies (B.7), it also satisfies (20). Thus, z is feasible to the ITND-2R model. We thus have that C˜2∗ ≤ C3∗ = C2∗ = C˜1∗ .  Appendix C. Supplementary data for Case II

Table C.5 The id, head, tail, mode, and length of each carrier link in the network for Case II. id

Head – Tail

Mode

Length

id

Head – Tail

Mode

Length

1 2 3 4 5 6 7 8 9 10

Harbin – Shenyang Shenyang – Dalian Shenyang – Beijing Beijing – Shijiazhuang Beijing – Tianjin Tianjin – Shanghai Dalian – Shanghai Shijiazhuang – Zhengzhou Yinchuan – Beijing Yinchuan – Lanzhou

Truck Truck Rail Rail Rail Sea Sea Rail Rail Rail

654.6 385.7 751.9 283.1 118.9 1229 1038 424.8 1251.9 487.3

25 26 27 28 29 30 31 32 33 34

Changsha – Chengdu Changsha – Zhuzhou Zhengzhou – Wuhan Wuhan – Zhuzhou Hangzhou – Zhuzhou Hangzhou – Shanghai Shanghai – Shenzhen Guiyang – Zhuzhou Guiyang – Kunming Kunming – Nanning

Rail Rail Rail Rail Rail Rail Sea Rail Rail Rail

1316.6 49.2 563.8 445.6 954 193 1531.7 808.5 572.2 909.1

(continued on next page)

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

103

Table C.5 (continued) id

Head – Tail

Mode

Length

id

Head – Tail

Mode

Length

11 12 13 14 15 16 17 18 19 20 21 22 23 24

Urumqi – Lanzhou Urumqi – Lhasa Xining – Lhasa Xining – Lanzhou Lhasa – Quxam Quxam – Chittagong Chittagong – Yangon Mandalay – Ruili Mandalay – Yangon Ruili – Kunming Lanzhou – Xian Chengdu – Xian Xian – Zhengzhou Chengdu – Kunming

Rail Rail Rail Rail Rail Rail Sea Rail Rail Rail Rail Rail Rail Rail

1948.1 2665.6 1904.4 198.6 707.3 1311.7 735.5 416.2 661.8 778.1 607.9 777.4 503.5 865.2

35 36 37 38 39 40 41 42 43 44 45 46 47 48

Nanning – Shenzhen Zhuzhou – Shenzhen Kunming – Shangyong Laocai – Kunming Hanoi – Laocai Hanoi – Langson Langson – Nanning Hanoi – Tay Ninh Tay Ninh – Bangkok Shangyong – Bangkok Yangon – Bangkok Bangkok – Singapore Shenzhen – Singapore Yangon – Singapore

Rail Rail Rail Rail Rail Rail Rail Rail Rail Rail Rail Sea Sea Sea

843.4 838.1 721 489.8 306.5 140.8 244.9 1391.5 805.5 1130 833 1557.1 2621.7 1234.3

Table C.6 O–D demands for Case II. Destination d

Origin o

1 5 7 10 12 16 17 20 27 30

14

23

29

31

32

34

16.98 23.63 2.52 0.00 0.00 35.66 0.00 0.00 0.00 0.00

30.12 41.96 4.51 14.43 3.31 63.13 40.22 13.53 9.24 1.09

0.00 0.00 0.00 2.41 0.55 0.00 6.64 2.26 1.56 0.19

0.00 0.00 0.00 0.29 0.07 0.00 0.00 0.00 0.00 0.02

0.00 0.00 0.00 1.64 0.37 0.00 0.00 0.00 0.00 0.00

143.28 199.12 21.53 69.44 15.89 299.20 190.82 64.90 44.89 5.39

References Alumur, S., Kara, B.Y., 2008. Network hub location problems: the state of the art. Eur. J. Oper. Res. 190 (1), 1–21. Alumur, S., Kara, B.Y., Oya, E.K., 2012. Multimodal hub location and hub network design. Omega 40 (6), 927–939. Ben-Akiva, M., Lerman, S.R., 1985. Discrete Choice Analysis: Theory and Application to Travel Demand. The MIT Press, Cambridge, MA. Bertsimas, D., Tsitsiklis, J.N., 1997. Introduction to Linear Optimization. Athena Scientific, Belmont, Massachusetts. Cascetta, E., 2009. Transportation Systems Analysis: Models and Applications, second ed. Springer, New York, NY. Crainic, T.G., Kim, K.H., 2007. Intermodal transportation. In: Barnhart, C., Laporte, G. (Eds.), Handbooks in Operations Research and Management Science, 14. North-Holland, Amsterdam, pp. 467–537. Croxton, K.L., Gendron, B., Magnanti, T., 2003. A comparison of mixed-integer programming models for nonconvex piecewise linear cost minimization problems. Manage Sci. 49 (9), 1268–1273. CRRC, 2012. Railway container freight billing method. http://www.crrcbj.com/oaGetData.do?fldId=53e4bd2a3b5eee42013b5fdf93c40 0 01. (in Chinese) (accessed: 14.11.15). Demir, E., Burgholzer, W., Hrušovský, M., Arıkan, E., Jammernegg, W., Van Woensel, T., 2016. A green intermodal service network design problem with travel time uncertainty. Transp.Res. Part B: Methodol. 93 (Part B), 789–807. Elhedhli, S., Hu, F.X., 2005. Hub-and-spoke network design with congestion. Comput. Oper. Res. 32, 1615–1632. Ernst, A.T., Krishnamoorthy, M., 1998. Exact and heuristic algorithms for the uncapacitated multiple allocation p-hub median problem. Eur. J. Oper. Res. 104 (1), 100–112. ESCAP, 2007. Toward an Asian Integrated Transport Network. Technical Report. United Nations. Kimms, A., 2006. Economies of scale in hub & spoke network design models: We have it all wrong. In: Morlock, M., Schwindt, C., Trautmann, N., Zimmermann, J. (Eds.), Perspectives on Operations Research: Essays in Honor of Klaus Neumann. Deutscher Universitätsverlag, Wiesbaden, Germany. LeBlanc, L.J., 1975. An algorithm for the discrete network design problem. Transp. Sci. 9 (3), 83–199. Liu, H., Wang, D.Z.W., 2015. Global optimization method for network design problem with stochastic user equilibrium. Transp. Res. Part B: Methodol. 72, 20–39. Loureiro, C.F.G., Balston, B., 1996. Investment selection model for multicontainer multimodal transportation networks. Transp. Res. Rec. 1522, 38–46. Luce, R.D., 1959. Individual Choice Behavior: A Theoretical Analysis. John Wiley & Sons, Inc., New York. Luna, H.P.L., Mahey, P., 20 0 0. Bounds for global optimization of capacity expansion and flow assignment problems. Oper. Res. Lett. 26 (5), 211–216. Meng, Q., Wang, X., 2011. Intermodal hub-and-spoke network design: incorporating multiple stakeholders and multi-type containers. Transp. Res. Part B: Methodol. 45 (4), 724–742. Meng, Q., Yang, H., Bell, M.G.H., 2001. An equivalent continuously differentiable model and a locally convergent algorithm for the continuous network design problem. Transp. Res. Part B: Methodol. 35 (1), 83–105. Meraklı, M., Yaman, H., 2016. Robust intermodal hub location under polyhedral demand uncertainty. Transp. Res. Part B: Methodol. 86, 66–85. Montesinos, V., Zizler, P., Zizler, V., 2015. An Introduction to Modern Analysis, first ed. Springer, New York, NY. Nocedal, J., Wright, S., 2006. Numerical Optimization, second ed. Springer, New York, NY. O’Kelly, M.E., 1987. A quadratic integer program for the location of interacting hub facilities. Eur. J. Oper. Res. 32 (3), 393–404. O’Kelly, M.E., Bryan, D., 1998. Hub location with flow economies of scale. Transp. Res. Part B: Methodol. 32 (8), 605–616. Serper, E.Z., Alumur, S., 2016. The design of capacitated intermodal hub networks with different vehicle types. Transp. Res. Part B: Methodol. 86, 51–65. Train, K.E., 2003. Discrete Choice Methods with Simulation. Cambridge University Press, Cambridge, UK.

104

X. Wang, Q. Meng / Transportation Research Part B 95 (2017) 76–104

Wang, S., Meng, Q., Yang, H., 2013. Global optimization methods for the discrete network design problem. Transp. Res. Part B: Methodol. 50, 42–60. Wang, X., Meng, Q., Miao, L., 2016. Delimiting port hinterlands based on intermodal network flows: model and algorithm. Transp. Res. Part E: Logist. Transp. Rev. 88, 32–51. Yamada, T., Russ, B.F., Castro, J., Taniguchi, E., 2009. Designing multimodal freight transportation networks: a heuristic approach and applications. Transp. Sci. 43 (2), 129–143.