Two-phase stochastic program for transit network design under demand uncertainty

Transportation Research Part B 84 (2016) 157–181

Contents lists available at ScienceDirect

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

Two-phase stochastic program for transit network design under demand uncertainty Kun An, Hong K. Lo∗ Department of Civil and Environmental Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong

a r t i c l e

i n f o

Article history: Received 19 September 2014 Revised 15 December 2015 Accepted 15 December 2015 Available online 13 January 2016 Keywords: Transit network design Service reliability Robustness Stochastic demand

a b s t r a c t This paper develops a reliability-based formulation for rapid transit network design under demand uncertainty. We use the notion of service reliability to confine the stochastic demand into a bounded uncertainty set that the rapid transit network is designed to cover. To evaluate the outcome of the service reliability chosen, flexible services are introduced to carry the demand overflow that exceeds the capacity of the rapid transit network such designed. A two-phase stochastic program is formulated, in which the transit line alignments and frequencies are determined in phase 1 for a specified level of service reliability; whereas in phase 2, flexible services are determined depending on the demand realization to capture the cost of demand overflow. Then the service reliability is optimized to minimize the combined rapid transit network cost obtained in phase 1, and the flexible services cost and passenger cost obtained in phase 2. The transit line alignments and passenger flows are studied under the principles of system optimal (SO) and user equilibrium (UE). We then develop a two-phase solution algorithm that combines the gradient method and neighborhood search and apply it to a series of networks. The results demonstrate the advantages of utilizing the two-phase formulation to determine the service reliability as compared with the traditional robust formulation that pre-specifies a robustness level. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Rapid transit systems, such as rapid bus, metro, light rail, are a lifeblood of cities. They enhance mobility and energy efficiency as well as mitigate roadway congestion. The Transit Network Design Problem (TNDP) is to decide the station locations, route alignments and frequencies of rapid transit lines (RTL) to serve the travel demands between specific origin-destination (OD) pairs (Fan and Machemehl, 2006). Transit services often run on fixed schedules which are published in advance. The schedules are not only crucial to achieve cost minimization in managing the fleet size, vehicle and crew scheduling, but also useful to commuters for departure time choices. Due to the stochastic nature of travel demand, it may happen that a heavy demand shows up on one day but a light demand on the other. Indeed, the optimal service network could be altered substantially by incorporating uncertain demand (Lium and Crainic, 2009). Developing a schedule that accommodates the uncertain demand variation in each day with a satisfactory level of comfort has thus attracted a lot of interest. Wang and Meng (2012) studied robust schedule design for liner shipping services considering the uncertainty in container handling time, demand, etc. Apart from the demand side, the uncertainty may arise from the supply side, such as link failure. De-Los-Santos et al. (2012) evaluated the robustness of rail transit network under two scenarios, with or without bridging ∗

Corresponding author. Tel: +852 2358 8389. E-mail address: [email protected] (H.K. Lo).

http://dx.doi.org/10.1016/j.trb.2015.12.009 0191-2615/© 2015 Elsevier Ltd. All rights reserved.

158

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

interruptions. Laporte et al. (2010) employed a game theoretic method for railway transit network design considering link failure and competing modes. In this paper, we consider the uncertainty arising from the demand side. The TNDP is generally formulated as a multi-objective optimization problem considering different objectives, such as (1) maximizing service coverage under budget constraint, (2) minimizing construction and operating costs while maintaining passenger travel time to a satisfactory level (Samanta and Jha, 2011). Route coverage and direct service provision are key elements in determining passengers’ mode choices (Hensher and Li, 2012); whereas construction and operating costs determine the economic feasibility or profitability of the system, and thus are major concerns to the operators. The current literature can thus be classified into two categories: one focuses on determining the rapid transit line alignment and station locations to maximize coverage and/or to minimize passenger cost (Yu et al., 2012; Laporte et al., 2005); the other aims at selecting the optimal station locations as well as frequencies with a fixed network topology to minimize the sum of user and operator costs (Wan and Lo, 2003, 2009; Laporte et al., 2007; Bruno et al., 1998). Specifically, Yu et al. (2012) developed an ant colony optimization model to maximize the demand density along a route, where transfers were specifically addressed as an important factor influencing passenger mode choices. Laporte et al. (2005) optimized the location of a single rapid transit line to maximize trip coverage. TNDP in the presence of competing modes, such as private car, is regarded as an extension of current studies. The interaction between the transit network to be constructed and the existing road network further adds complexity to the problem. Wan and Lo (2003, 2009) formulated the multi-modal network design problem. Recently, Perea et al. (2014) proposed a mixed integer non-linear programming model to add a new railway station and a new road junction on a road-rail network in the presence of modal competition. While the deterministic TNDP has been explored extensively, few studies have investigated the issue of demand uncertainty associated with this problem. Typically there are two methods to address stochastic demand: one is stochastic optimization which aims at minimizing the expected cost by assuming that the demand follows a certain probability distribution; the other is robust optimization which aims at minimizing the cost associated with the worst case scenario (An and Lo, 2015). Robust optimization requires less information on the stochastic demand than stochastic optimization. It assumes that the demand distribution is only partially known, such as first and second moments, and belongs to a given family of probability distribution. One may refer to Gabrel et al. (2014) for a review of recent developments in robust optimization. Stochastic programming typically draws upon a certain sampling method and adopts the sample average to approximate the cost expectation, which can be expressed as a large-size mixed integer linear programming (MILP). It can be solved by commercial software such as CPLEX or by the L-shaped method (Slyke and Wets, 1969; Birge and Louveaux, 1988). Owning to the formidable computational work by exact methods, many studies resort to heuristic algorithms for solving large-size problems, such as neighborhood search methods, genetic algorithms and hybrid search methods (Guihaire and Hao, 2008). Robust optimization focuses on the worst case scenario, and thus is of a similar level of complexity as its deterministic counterpart. The downside is that its solutions may be conservative. Some studies thus turn to redefining the uncertainty set such that all possible realizations within the set are considered, whereas those outside are ignored. In the end, it is important to decide the size of the uncertainty set over which the worst scenario is generated, which ideally should strike a proper balance between the resultant system cost and the level of robustness achieved (Ben-Tal and Nemirovski, 1999; Bertsimas and Sim, 2004). As the size of the uncertainty set is decisive to the solution quality, choosing it inappropriately may, on one hand, fall short of the ability to hedge against uncertainty or, on the other, produce an overly conservative solution. Developing a method to rationally determine the size of the uncertainty set, which would relieve the dependency on expert opinion, is non-trivial. In response to this question, some studies use the cost variance across scenarios as a measure of solution robustness and incorporate it into the objective function. This approach allows an explicit mean-variance tradeoff. List et al. (2003) examined the bus fleet sizing problem via a robust optimization model by incorporating the meanvariance tradeoff in the objective function under demand and vehicle productivity uncertainty. Yan et al. (2012) developed a robust optimization model for scheduling of a fixed bus route aiming at minimizing the sum of the random schedule deviations. These studies typically focus on bus travel time reliability. Both the mean and variance of travel time are important attributes to measure system performance; hence the objective of minimizing a combination of the mean and variance of travel time is sensible. In this study, however, we aim at providing transit services to address demand uncertainty; cost variance, while interesting to know, does not quite capture the objective of the operator to minimize the expected cost. Robust and stochastic optimizations are two complementary ways to deal with uncertainty, each having its own advantages and disadvantages. In spite of their differences, researchers have made renewed efforts to bridge these two approaches. Bertsimas and Goyal (2010) quantified the performance of robust optimization with different distributions. Chen et al. (2007) investigated stochastic programming from the perspective of robust optimization. In this paper, we provide a stochastic perspective on robust optimization. Namely, instead of fixing the size of the uncertainty set a priori, when more information about the underlying distribution is available, can we limit the size of the uncertainty set in a rational manner so as to strike a proper balance between robustness or reliability of the system and its resultant cost? Such an analysis will be helpful in identifying the value of collecting information about the underlying probability distribution. In this paper, we develop another way to evaluate the outcome of the uncertainty set in TNDP. Recently, we formulated the TNDP under demand uncertainty for system optimal flows by considering the combination of two services types, (i) rapid transit lines (RTL) or regular services and (ii) demand responsive or flexible services. Defining the notion of service reliability (SR), we proposed a two-phase model to separate the otherwise intertwined decisions over the deployment of these two service types (Lo et al., 2013; An and Lo, 2014a, 2014b). The RTL are designed to cover the stochastic demand up to a certain specified SR. The RTL operate on dedicated right-of-ways or tracks with fixed schedules, whereas flexible services

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

159

Fig. 1. Illustration of the station gate sub nodes, line sub nodes, in-station links, out-station links, dummy origins, dummy destinations, and dummy links. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

are subcontracted to a third party to carry the demand overflow on a particular day. The transit operator jointly plans the alignments and schedules of the RTL and the extent of deploying third party flexible services. To service operators, flexible services have higher per unit costs; hence, their objective is to determine the optimal combination of these two service types to minimize the overall construction and expected operating cost. The flexible services cost can also be viewed as the penalty of demand loss if only regular transit services are provided. The SR chosen in designing the regular services implies a certain level of expected cost in deploying the flexible services. In this way, the SR is not arbitrarily set but is linked to explicit cost terms in deploying the flexible services. And by optimizing the overall combined cost of the two service types, the associated optimal SR is also internally determined. In this study, we extend the consideration to both system optimal and user equilibrium passenger flows, and develop a two-phase solution algorithm that combines the improved gradient method and neighborhood search. This paper is structured as follows: Section 2 develops the model formulations and solution algorithms for system optimal (SO) and user equilibrium (UE) passenger flows. Section 3 applies the mathematical formulations and algorithms to a series of networks, and compares the solutions obtained by the two-phase model and the conventional robust model. Finally, Section 4 concludes the paper with some possible future research directions. 2. Model formulation 2.1. Problem setting We consider a given transportation network G(N, A) defined by a node set N and a link set A. Stations i and j, i, j ∈ N, i  j are connected by link (i, j) ∈ A with distance lij . This paper aims at determining the RTL alignments Y and frequencies = f for up to Rmax lines, as well as the deployment of flexible services Z at the same time in a most cost efficient way while satisfying the stochastic demand. R = {r} denotes the set of RTL. fmin , fmax , respectively, denote the lower and upper bound of line frequency. We assume each transit line is acyclic but provides services for two directions with the same frequency. Recently, Walteros et al. (2015) decomposed stations into gate sub-nodes and route sub-nodes to represent feasible routes. Similarly, to represent passenger transfer cost, alighting and boarding within stations, we decompose each station i into sub-nodes, denoted by ir :r = 0 for station gate at node i (to capture demands arriving or departing at the station) and r = 1, 2 . . . Rmax for line sub-nodes at station i (to capture demand transfers between lines at station i). In-station links (iu , iv )∀i ∈ N, u, v ∈ R∪{0} and u = v are constructed to fully connect the station and line sub-nodes within a station, with its flow representing (1) the amount of passenger transfers between lines if u, v ∈ R, and u = v; (2) the amount of passengers boarding to line v if u = 0, v ∈ R; (3) the amount of passenger alighting from line u if u ∈ R, v = 0. Out-station links (iu , ju ) ∀(i, j) ∈ A, u ∈ R are constructed to connect line sub-nodes in different stations with the same line index, with its flow representing passengers traveling from station i to j on line u. Note that out-station links are only constructed for feasible links (i, j) ∈ A. To select the origin and destination station for each transit line flexibly, we introduce a dummy origin set S = {Sr , r ∈ R} and dummy destination set T = {Tr , r ∈ R}, so that every line r has a fixed dummy origin Sr and dummy destination Tr . Dummy links (Sr , ir ), i ∈ N, r ∈ R and (ir , Tr ), ∀i ∈ N, r ∈ R are constructed to connect every line sub-node ir , ∀i ∈ N, r ∈ R with the dummy origin and destination, respectively. A graphical illustration for a three-node network is given in Fig. 1 to show the related station sub-nodes, line sub-nodes, in-station links, out-station links, dummy origin and destination nodes and dummy links. Three location nodes A, B and C in the network, as represented by the blue ellipses in Fig. 1, are considered to be candidate stations for up to two transit lines. The subscripts represent the sub-node index. Solid lines

160

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

without arrows indicate that both directions are available. One dummy origin and one dummy destination are constructed for each transit line, as represented by the dark grey nodes. Dummy links are represented by dotted lines with arrows indicating their directions to connect to the dummy nodes and each line sub-nodes. By treating the nodes Sr , r = 1 or 2 as origins and nodes Tr , r = 1 or 2 as destinations, it requires the transit lines to traverse through the candidate stations in the actual network. The node connected to a Sr node will become the starting station of the transit line; similarly, the one connected to a Tr node will become the terminal station of that transit line. The horizontal lines connecting A1 ,B1 ,C1 or A2 ,B2 ,C2 are the out-station links for line 1 or 2, representing the possible transit line segment. Take station A for example, it contains one station sub-node, represented by the light grey node A0 , and two line sub-nodes A1 and A2 for RTL lines 1 and 2, respectively. The links connecting A0 with A1 and A2 enable passengers to board or alight from the RTL, and the link connecting A1 and A2 enables passengers to transfer between the RTL lines 1 and 2. The demands are loaded to the network through the station sub-nodes A0 ,B0 ,C0 if a station is the origin of an OD pair d ∈ D. Passengers from station A to C can take a route A0 -A1 -B1 -B2 -C2 -C0 if link A to B is covered by line 1 and link B to C is covered by line 2. Passengers have to transfer from line 1 to line 2 at station B. For ease of referencing, all the variable definitions are provided in Appendix. The assumptions are as follows: (a) (b) (c) (d) (e)

The probability distributions of the stochastic demands can be estimated. Flexible services, once deployed, provide door to door services for OD pairs d ∈ D. Road congestion of flexible services is not substantial. Passengers do not transfer between transit lines and flexible services. A link can be used by multiple transit lines.

It is common to consider that passenger arrivals follow a certain probability distribution. For an existing line, the passenger arrival distribution can be obtained by regression analysis of historical data. For a newly developed line, it can be obtained by drawing upon data of other cities with similar land use pattern and population. (b)–(d) define the function and operation mechanism of the flexible services, which are initiated by the transit operator and operated by a third-party to serve the demand overflow. Flexible services run on the shortest paths between OD pairs on the roadway network. Due to their shorter travel times and zero waiting time, flexible services are passengers’ first choices once provided. To the transit operator, these two services have different costs. c) assumes that flexible services should not significantly aggravate the roadway congestion. However, the higher unit cost of flexible services to the company somewhat discourages its widespread use. e) describes that rail track can be shared by multiple lines. The decision variables include the following. A set of binary variables Y = {Yir j } that denotes whether an out-station link u v is on a transit line. Since out-station links (iu , jv ) are defined in coordination with line r, i.e., u = v = r, the subscripts u, v in (iu , jv ) are thus redundant and omitted for transit line decision variables Yir j in the following description. Yirj is 1 if link (ir , u v jr ) is on line r; 0 otherwise. f = { fr } denotes the set of frequency variables. To formulate a linear program, we define a new variable set y = {yrij } whose elements are defined as the product of Yirj and fr , i.e. yrij = Yirj ∗ fr . The binary variable Y i j equals

1 if link ij is covered by any line r ∈ R, zero otherwise, with its set denoted by Y = {Y i j }. The binary variable set W = {Wi } determines whether a station i is on a transit line: Wi = 1 if any transit line r traverses station i, 0 otherwise. Station links are constructed automatically if Wi = 1. Note that Y, y, W can be uniquely determined by Y and f.Xid j represents the passenger u

v

flow on RTL for OD pair d with its set denoted by X = {Xid j }. The subscripts (iu , jv ) categorize Xid j into four subsets: (1) u

v

u

v

d for i = j, u = 0, v ∈ R, which represents the passenger flow boarding to transit line v at Boarding passenger flow subset Xboard d for i = j, u ∈ R, v = 0, which represents the passenger flow alighting from station i; (2) Alighting passenger flow subset Xalight

d for i = j, u, v ∈ R, u = v, which represents the amount transit line u at station i; (3) Transfer passenger flow subset Xtrans f er d of passengers transferring from line u to v at station i; (4) In-vehicle passenger flow subset Xin− vehicle for i = j, u = v ∈ R,

which represents the passenger flow traveling from station i to j on line u. Zd is the amount of flexible services provided for OD pair d and is also equal to the amount of passengers taking this service, with its set denoted as Z = {Z d }. Note that the decision variables are directional. The problem is to decide for each transit line its forward direction, then services for the backward direction is added automatically. In the following, we first state the robust formulation for the TNDP with a predetermined uncertainty set, then we describe the two-phase reliability-based formulation for SO and UE passenger flows, respectively. 2.2. Robust formulation with a predefined uncertainty set 2.2.1. Robust formulation with system optimal (SO) flows In the domain of robust optimization, the approach typically seeks to search for the best solution that is feasible to all the possible realizations of the random parameters. Robust TNDP can be regarded as the optimization of the system cost under the worst case scenario. A broader uncertainty set enhances the level of robustness, which would require a higher RTL capacity to hedge against demand fluctuation, thus is beneficial to passengers, but detrimental to the operator’s benefits. Let Bd be the random demand with mean bd and standard deviation σ d , and its cumulative distribution function (CDF) denoted by  d . If the random demand follows a continuous distribution  d , defined from zero to infinity, and all possible demand realizations within  d were to be met, then this would require a service with an infinite capacity. To overcome this

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

161

problem of providing overly conservative solutions, typically the approach utilizes a redefined uncertainty set that excludes the events with low probabilities. Let d = [bd − n∗ σ d , bd + n∗ σ d ] be the redefined uncertainty set where n determines its breadth. Denote N = N ∪ S ∪ T as the set of all nodes including the dummy origin and destination nodes and A as the set of all links including the dummy links. The traditional robust formulation under SO is presented as follows:

(P0 − SO )

min

max

f,Y,y, Y,W,Z B∈,X(B )

φ=



ci1j yrij +

r∈R i j∈A

+





ci2jY i j +

i j∈A

d 5 ci4j Xin− vehicle +c

i j∈A d∈D



c3Wi

i∈N

  1  1 d τ+ Xboard + c5 Xd 2 fr 2 fr trans f er r∈R d∈D

(1)

r∈R d∈D

Subject to



Yirj ≤ 1,

∀r ∈ R, i ∈ N ∪ S

(2)

Y jir ≤ 1,

∀r ∈ R, i ∈ N ∪ T

(3)

j∈N

 j∈N



Yikr =

i∈N



Ykrj ,

∀r ∈ R, k ∈ N

(4)

j∈N

Yirj + Y jir ≤ 1, ∀i j ∈ A, r ∈ R   Yij ≤ Yirj + Y jir , ∀i j ∈ A and Yij ≥

r∈R Yirj + Y jir ,

(5)

i< j

∀r ∈ R, i j ∈ A and i < j ∀i ∈ N, j ∈ N

Wi ≥ Y i j , Wi ≥ Y ji ,

(6) (7) (8)

Yirj = 0

or

1,

∀r ∈ R, i j ∈ A

Yij = 0

or

1,

∀i j ∈ A

(10)

Wi = 0

or

1,

∀i ∈ N

(11)

fmin ≤ fr ≤ fmax , ∀r ∈ R   fmin ≤ yrij + yrji ≤ fmax ,

(9)

(12)

∀i j ∈ A

(13)

r∈R

∀i j ∈ A, r ∈ R

yrij − fr ≤ 0, yrij − MYirj ≤ 0,





(14)

∀i j ∈ A, r ∈ R

M Yirj − 1 − yrij + fr ≤ 0,

(15)

∀i j ∈ A, r ∈ R

yrij = frYirj , ∀i j ∈ A, r ∈ R  Yirj ≤ || − 1, ∀ ∈ N, || ≥ 2, r ∈ R i∈ j∈

Xid0 ir −

u∈R∪{0}







(17) (18)



Bd , if i is the origin of OD d −Bd , if i is the destination of OD d , ∀i ∈ N, d ∈ D r∈R r∈R 0, otherwise     d d Xir iu − Xiu ir + Xidr jr − X jdr ir = 0, ∀i ∈ N, r ∈ R, d ∈ D



(16)

Xidr i0 =

u∈R∪{0}



Xidr jr ≤ C yrij + yrji ,

j∈N

(19) (20)

j∈N

∀i j ∈ A, r ∈ R

(21)

d∈D

ci1j , ci2j , ci4j represent the transit link operating cost, transit link construction cost, passenger travel cost on link ij, respectively, which are all proportional to the link distance lij . c3 denotes the station construction cost. c5 represents the passenger value of time. To make the notations succinct, we omit the subscript of passenger flow variable Xid j in the objective function u v but use its category name instead. The objective function is comprised of six cost components, i.e. transit operating cost, transit line construction cost, station construction cost, passenger in-vehicle travel time cost, passenger waiting time cost

162

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

for boarding, and passenger transfer penalty. The passenger waiting time to get onboard is set as half of the headway of the arriving transit line (Yu et al., 2012), as shown by the fourth expression. The last expression calculates the transfer cost consisting of a transfer penalty τ and waiting time cost. Constraints (2) and (3) ensure that only one station is directly connected to node i from upstream and downstream, respectively, if node i is on line r. Constraint (2) also indicates that at most one link can be generated from the dummy origin Sr . Constraint (3) ensures at most one link is ended at destination Tr . (4) is the connectivity constraint to ensure that there are exactly two links connecting each station on line r. (5) represents that at most one direction of a link can be occupied by one transit line. Constraints (6) and (7) ensure that a rail track is constructed on link ij if any line r ∈ R goes through it. Constraint (8) ensures that a station is constructed if any line traverses through it. Constraint (13) states that the frequency of link segment ij is restricted by the maximum and minimum frequency boundaries. Constraints (14)–(16) are linear representations of the nonlinear Eq. (17). (17) is shown here merely for illustration purposes and will be removed in the solution process. M is a very large positive number. The cycle breaking Constraint (18) states that the number of links connecting nodes in the node sub set ,  ∈ N, must be strictly less than the node number. In a cyclic route, the number of links and nodes is equivalent. This constraint is not added unless some cycle solution is found. Constraints (19) and (20) represent the passenger flow balancing conditions for station sub-nodes and line sub-nodes, respectively. The first (second) summation in Constraint (19) calculates the amount of passengers boarding (alighting) from station node i0 . Constraint (19) states that the passenger flow entering into (leaving) the transit system through station i is equal to the demand if i is the origin (destination) of one OD pair. Similarly, Constraint (20) represents that the amount of passengers on the incoming links of a line sub-node is equal to passenger flow on its outgoing links. Even though the right hand side of (21) sums the forward and backward directions of rapid transit link ij frequencies, i.e. yrij + yrji , it can only contribute to the constraint once. The reason is that if Yirj is one, then Y jir must be zero for acyclic lines. The nonlinear term 1/fr involved in the passenger transfer and boarding costs makes P0-SO a non-convex non-linear mixed integer programming with uncertainty constraints, which is very difficult to solve. Hence in the case study, we select fr = ( fmin + fmax )/2 as an initial point to replace fr in passenger waiting and transfer cost terms in the objective function. P0-SO is thus relaxed into a min–max problem with linear constraints and objective function. After solving this relaxed P0-SO, we get the optimal frequency fr for this relaxed problem. We then update fr := fr + 12 ( f r − f r ) and repeat the previous steps until | f r − f r | ≤ ε for all r ∈ R. The final optimal frequency is taken as 12 ( f r + f r ). This relaxed P0-SO is a min–max problem in which the random demand is confined by the uncertainty set B ∈ . Given a design (f, Y, y, Y, W, Z ), it searches for the worst-case scenario B ∈  which has a maximum total system cost. Then it seeks to minimize this cost by modifying the design. Different linearization methods have been proposed to tackle this problem for various forms of the uncertainty set (Ben-Tal and Nemirovski, 1999; Bertsimas and Sim, 2004). In P0-SO, we can relax the equality constraints (19) by these inequality constraints:



Xid0 ir −

r∈R



Xidr i0

r∈R

⎧ d ⎨≥ B ,

if i is the origin of OD d if i is the destination of OD d , otherwise

≤ −Bd , ⎩ = 0,

∀i ∈ N, d ∈ D, B ∈ 

(22)

without affecting the optimal solutions, since a larger demand shall incur an equal or larger cost under the principle of system optimal flows (Ben-Tal et al., 2011). In (22), the random demand is confined by an uncertainty set B ∈ . This constitutes a robust counterpart that seeks to search for an optimal solution (f, Y, y, Y, W, Z, X ) that is feasible for any demand realization in . The robust counterpart can thus be written as:

(P0 − SO’ ) min

,W,Z f,Y,y,Y

+ c5

φ=



ci1j yrij +

r∈R i j∈A

 r∈R d∈D

τ





ci2jY i j +

i j∈A



c3Wi +

i∈N



d 5 ci4j Xin− vehicle + c

 1 Xd 2 fr board r∈R d∈D

i j∈A d∈D

1 + Xd 2 fr trans f er

s.t. (2 ) − (18 ), (20 ) − (21 ),

B − AX ≤ 0,

∀Bd ∈ [bd − n∗ σ d , bd + n∗ σ d ],

(23)

where A is the parameter matrix on the left hand side of (22). Since the uncertainty only appears in the right hand side vector B in (22), it is a special case of the robust formulation in Soyster (1973). Employing the same technique in Bertsimas and Sim (2004), the stochastic demand must take its largest value, i.e. Bd = bd + n∗ σ d , such that B − AX ≤ 0 holds for any Bd ∈ [bd − n∗ σ d , bd + n∗ σ d ]. That is, due to the special structure of (22), for independent stochastic demand bounded by a box uncertainty set, one can substitute the random demand by its upper bound. 2.2.2. Robust formulation with user equilibrium (UE) flows

(P0 − UE )

min

max

f,Y,y, Y,W B∈,X(B )

 r∈R i j∈A

ci1j yrij +

 i j∈A

ci2jY i j +

 i∈N





c3Wi + t f, Y, y, Y, W + g(f, Y, y, Y, W )

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

163

Subject to (2)–(18) and X is solution of the lower level problem:

min t (f, Y, y, Y, W ) Subject to (19 ) − (21 ), X

where

t (f, Y, y, Y, W ) =



d 5 ci4j Xin− vehicle +c

r∈R d∈D

i j∈A d∈D

and

g(f, Y, y, Y, W ) =



  1  1 d τ+ Xboard + c5 Xd 2 fr 2 fr trans f er

−c5 βirj

i j∈A r∈R



r∈R d∈D

Xidr jr

d∈D

P0-UE is a bi-level min–max problem where the UE passenger flows are calculated in the lower level problem. To include the effect of the hard capacity constraints, we introduce an overflow delay for transit line segments operating at capacity. This overflow delay includes passenger discomfort on a train fully packed to its capacity. Consequently, passenger’s cost has two components: free-flow path travel time cost and overflow delay if a RTL line segment is operated at capacity. We take the negative Lagrange multiplier −βirj of (21) as the overflow delay of RTL line segment operating at capacity (Meng et al., 2008; Lam et al., 1999). Through minimizing the passenger travel time cost t (f, Y, y, Y, W ), i.e. solving the linear program of the lower level problem generates the UE passenger flows (An and Lo, 2014b; Lam et al., 1999). Note that paradox could occur under UE in the sense that a smaller demand could incur a larger system cost. In order words, the maximum demand combination may not constitute the worst case scenario, and the set of equality Constraint (19) cannot be relaxed by the set of inequality Constraint (22). This imposes great computational difficulty in searching for the actual worst case demand scenario under UE with hard capacity constraints. One possible way is to extend the cutting plane method (Lou et al., 2009; Chung et al., 2012) designed for robust discrete network design problem without hard capacity constraints to our model. Nevertheless, the cutting plane method involves a tedious process of cutting and resolving the formulation for different demand combinations. Moreover, the hard capacity constraints needed for this problem add further complexity. To the best of our knowledge, developing an efficient solution procedure for the problem of robust network design under both UE and hard capacity constraints is yet to be tackled in the literature. It is a worthwhile endeavor and left for our future research. As such, we merely formulate the min–max bi-level problem here for the sake of completeness. In this study, we propose an alternative formulation of the problem through the notion of service reliability and solve it accordingly, as will be discussed in Section 2.3. The solution quality under both SO and UE depends on the selection of n, i.e. the size of the uncertainty set. It hence raises the question of choosing an appropriate n to trade off the level of robustness and the resultant system cost. In the following sections, we make use of the service reliability (SR) ρ as a new variable to separate the formulation P0-SO into two phases. Note that ρ is a vector, whose elements may be different for different OD pairs. ρ works as an internal variable to link the two phases. In phase-1, given a certain specified ρ, we determine the RTL alignments Y and frequencies f to meet the stochastic demand up to a certain service reliability ρ. To capture the cost incurred by demand realizations beyond ρ, we propose flexible services to carry the demand overflow that exceeds the capacity of RTL. The flexible services deployment, which may vary from day to day, is determined in the phase-2 problem. Qualitatively speaking, a high SR indicates a RTL system with a high capacity to accommodate demand fluctuations. The downside is that when the realized demand is low, the RTL capacity may not be fully utilized, implying a certain level of waste. In contrast, a RTL system designed with a low SR needs to activate flexible services more frequently with higher costs. The SR is, thus, a variable that trades off the RTL and flexible services provision to minimize the expected total cost. In this formulation, passenger assignment on these two services is realized in the second phase, which generates two schemes: one is system optimal assignment, i.e. passengers are assigned to routes to minimize the total system cost; the other is user equilibrium, i.e. all RTL passengers on the same OD pair are assigned to routes with the same minimal cost. We study the influences of passenger assignment patterns on the transit network configurations separately in Sections 2.3 and 2.4. 2.3. Stochastic formulation with system optimal flows 2.3.1. Model formulation In the stochastic formulation, with the estimated demand distributions, the primary objective is to find the optimal SR that minimizes the expected total system cost. Passengers’ preferences over the transit and flexible services are not considered in the SO passenger assignment. In phase 1, we determine the location of RTL stations and frequencies that can cover the demand up to B¯ d = d−1 (ρ d ), ∀d, where  d is the CDF of the random demand. The passenger in-vehicle traveling time cost on RTL is included in the objective function. Doing so will guide the model for a RTL solution that considers both passenger and operator costs. (P1) Phase-1:

min

f,Y,y, Y,W

 r∈R i j∈A

ci1j yrij +

 i j∈A

ci2jY i j +

 i∈N

c3Wi +

 i j∈A d∈D

d ci4j Xin− vehicle

(24)

164

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

Subject to (2)–(16), (18) and (20)–(21)



Xid0 ir −

r∈R



⎧ d ⎨B¯ ,

Xidr i0 =

r∈R

B¯ d = d−1 (ρ d ),

if i is the origin of OD d if i is the destination of OD d , otherwise

−B¯ d , ⎩ 0,

∀i ∈ N, d ∈ D

d∈D

(25)

(26)

In the passenger flow conservation Constraint (25), the demand to be covered by RTL B¯ d is fixed by (26), whereas in P0-SO’, it is a random parameter Bd constrained by an uncertainty set d , i.e. Bd ∈ d . ρ d , to some extent, determines the upper bound of d and hence assists in transforming the stochastic problem P0-SO’ to a deterministic problem P1 with fixed demand. The transformation alleviates the cumbersome linearization and computation. After solving P1, the transit cost is recalculated as:



J (ρ ) =

ci1j yrij +

r∈R i j∈A



ci2jY i j +

i j∈A



c3Wi

(27)

i∈N

The passenger cost in (24) is an estimation by fixing the demand at B¯ d , and is abandoned in (27) to avoid double counting of the passenger cost which will be more accurately determined in phase-2. The actual passenger and flexible services cost subject to demand realizations are evaluated in the phase-2 problem, which aims to minimize the expected cost of flexible services and passenger cost. To address the stochastic demand, a sampling procedure is adopted to estimate the expected cost terms by their sample averages. One may choose the appropriate sampling method as well as the sample size according to the specified problem. Assume that m scenarios with equal probability are generated, each scenario e ∈ E with 1 . The subscript or superscript e for all variables below denotes the scenario e. With the RTL alignments a probability of pe = m Y and frequencies f fixed by solving P1, phase-2 determines the deployment of flexible services. (P2) Phase-2:

min Qe = X,Z



cd6 Zed +

Xid,e 0 ir

r∈R



⎧ d ⎨Be − Zed ,

Xid,e r i0

=

u∈R∪{0}

Xid,e − r iu

Zed

⎩

−

Bde ,

0,

 u∈R∪{0}



Xid,e + u ir

 j∈N



Xid,e ≤ C Yirj + Y jir fr , r jr

  1  1 d,e 5 X + c + X d,e ∀e ∈ E τ 2∗ fr board 2∗ fr trans f er r∈R d∈D

i j∈A d∈D

r∈R





−

d,e 5 ci4j Xin− vehicle +c

d∈D

Subject to





(28)

r∈R d∈D

if i is the origin of OD d if i is the destination of OD d , otherwise

Xid,e − r jr



X jd,e =0 r ir

∀i ∈ N, d ∈ D

∀i ∈ N, r ∈ R, d ∈ D

(29)

(30)

j∈N

∀i j ∈ A, r ∈ R

(31)

d∈D

cd6 denotes the combined flexible services cost for OD pair d consisting of the company operating cost and passenger travel cost. In the case study, for results exposition, we will separately report these two cost components. With Y and f predetermined in phase-1, P2 becomes linear program problem that can be solved much more efficiently. The passenger flow can be carried by RTL are restricted by its capacity y = {yrij = Yirj ∗ fr } as reflected in (31). The additional capacity provided by the flexible services Zed , ∀d, tailored for the demand overflow beyond ρ, is calculated for each scenario e ∈ E in phase-2. The expected phase-2 cost is calculated as the weighted sum of all scenarios as shown in (32).

Q (ρ ) =



pe Qe

(32)

e∈E

We then combine the deliberately separated two phases to formulate the overall system cost:

min

X,Y,Z,f,ρ

φ (ρ ) = J (ρ ) + Q¯ (ρ )

(33)

For a specified SR ρ, the RTL cost J(ρ) in (33) is solved and fixed in phase-1, whereas the expected flexible costs and the passenger cost Q¯ (ρ ) in (33) are solved in phase-2. By optimizing ρ, we minimize the combined system cost. We develop a two-phase solution procedure over ρ to solve the phase-1 and phase-2 problems sequentially, as will be described below. We should note that the combined phase-1 and phase-2 problem is non-convex, as are typical for two-stage stochastic programs. Hence, the optimal ρ such determined constitutes a local, rather than global, optimum.

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

165

2.3.2. Two-phase solution procedure The proposed two-phase solution procedure combines a neighborhood search with the gradient procedure proposed by Lo et al. (2013). In addition, the two-phase method here utilizes an analytical approach to obtain the gradient instead of the original trial-and-error method, which greatly improves the efficiency of the existing gradient method. The algorithm starts from an arbitrary initial ρ, and searches for its decent direction iteratively for the optimal ρ∗ so as to minimize the expected overall system cost (33). In our recent studies, a perturbation analysis was adopted by Lo et al (2013) and An and Lo (2014b) to calculate the partial derivative of the total system cost with respect to ρ. The main characteristics are summarized as follows. Firstly, we change one element in ρ = (ρ 1 , ρ 2 , . . . ρ d , . . . , ρ dmax ) by a small positive number δ d to obtain a trial value ρd = (ρ 1 , ρ 2 , . . . , ρ d + δ d , . . . , ρ dmax ). dmax is the number of OD pairs, also the dimension of ρ. Then we solve the phase-1 φ (ρd )−φ (ρ )

is utilized as the and phase-2 problems for the trial ρd and obtain the total system cost φ (ρd ). The cost ratio δd partial derivative. We then repeat the whole process for each element of ρ to formulate the gradient of system cost φ over ρ. As the phase-1 and phase-2 problems are solved repeatedly for every element of ρ, the computation time highly depends on the dimension of ρ, which is equal to the number of OD pairs. To improve the computational efficiency, in this study, we take advantage of the dual solutions of the phase-2 problem to obviate the tedious numerical calculation of the previous approach. Specifically, during the gradient calculation in our previous study (An and Lo, 2014a), an interesting phenomenon is observed; i.e. when the deviation δ d in ρ was small enough (say, less than 0.1), the RTL alignments Y typically remain unchanged. Only the frequencies f are modified to accommodate the demand change in P1, since RTL realignment will result in substantial system cost changes whereas changes in frequencies do not. Thus a miniscule perturbation in demand typically would not influence the RTL alignments. Inspired by this observation, we assume that a small change in ρ d only modifies the operating frequencies f, but not the RTL alignments Y. Hence, the original P1 is relaxed to a simple LP, which substantially reduces the computation time. Furthermore, since f only appears in the right hand side of (31), we can simply use the dual solution of P2 to obtain the cost difference, which emancipates the computational intensive sampling effort of the phase-2 problem. Proposition 1 below details the calculation. Proposition 1. Suppose (ρ, Y, f, β , J (ρ ), Q (ρ ), φ (ρ )) is the optimal solution of the phase-1 problem and phase-2 problem  } is the set of dual solution expectations, and the element βir,e is the for a specified SR ρ, where β = {βirj : βirj = e∈E pe βir,e j j

dual solution corresponding to Constraint (31) in the phase-2 problem for link ij on line r of scenario e. (ρd , Y, fd , J(ρd )) denotes the optimal solution of phase-1 after a small perturbation δ d in ρ while maintaining the same RTL alignment Y. The frequency difference set is defined as fd = fd − f = { frd : frd = frd − fr }. Then the phase-2 cost difference can be calculated by:

Q (ρ d ) − Q (ρ ) =







C Yirj + Y jir frd βirj

(34)

i j∈A r∈R

The system cost difference is:

  φ (ρ d ) − φ (ρ ) = J ρ d − J (ρ ) + Q (ρ d ) − Q (ρ )

(35)

And the gradient of system cost with respect to SR is:

φ ( ρ ) φ ( ρ d ) − φ ( ρ ) = ρ d δd

(36)

Proof. A sufficiently small perturbation δ d in ρ induces a small change frd in frequency. With the assumption that the small perturbation will not change the line alignments, hence C (Yirj + Y jir ) on the right hand side of (31) maintains unchanged. For a specific scenario e, making use of the shadow price, the phase-2 cost difference is the dual solution mul  r r d r,e tiplied by the change in right hand side value, i.e. i j∈A r∈R C (Yi j + Y ji ) f r βi j . Combining all the scenarios generated,    . And we get Q (ρd ) − Q (ρ ) = the expected phase-2 cost difference is Q (ρd ) − Q (ρ ) = e∈E pe i j∈A r∈R C (Yirj + Y jir ) frd βir,e j    r,e r r d i j∈A r∈R C (Yi j + Y ji ) f r e∈E pe βi j by moving pe inside. The system cost difference and the gradient hence can be obtained accordingly.  To summarize, we obtain the dual solutions associated with Constraint (31) for each scenario and calculate the dual  first. The phase-2 cost difference required in the gradient calculation can be obtained by (34) instead expectation e∈E pe βir,e j of resorting to repeatedly solving a large number of LPs in the phase-2 problem. The computation time is thus reduced to 1/ (1+dmax ), where dmax is the number of OD pairs and the dimension of ρ as compared with the gradient method developed in our earlier studies. We provide the searching procedure for the optimal ρ∗ in the following table. Let k be the iteration number, kmax be the maximum iteration number. The superscript ∗ of each variable denotes the optimal solution.

166

1 2 3 4 5 6 7 8 9 10 11 12 13

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

Initializeρk = (ρkd ), ∀d, k = 1. For (k = 1; k < kmax ; k + +) { Solve the phase-1 problem P1 to obtain the rapid transit alignment and the corresponding frequencies (Yk , fk ) for the given ρk . For (e=1 to m) {Solve the phase-2 problem P2 to obtain the flexible services deployment Ze,k and calculate the dual solution expectation β k }. Calculate the expected phase-2 cost Q (ρk ) through (32). Calculate the overall system cost φ (ρk ) through (33). If (φ (ρk ) < φ ∗ ), update the optimal solution (φ ∗ , ρ∗ , Y∗ , f∗ ). If (|φ (ρk ) − φ (ρk−1 )| < ε ), Break, take the optimal solution as (φ ∗ , ρ∗ , Y∗ , f∗ ). Determine ρ for the next step. Denote the solution set obtained in current step k as (Yk , fk , β k , φ (ρk )) For (d=1 to dmax ) { Set ρkd = (ρ 1 , ρ 2 , . . . , ρ d + δ d , . . . , ρ dmax ) Solve phase-1 problem P1 to obtain frequency difference fd with given ρkd and Yk . Calculate phase-2 cost difference Q (ρkd ) − Q (ρk ) through (34). Calculate the system cost difference φ (ρkd ) − φ (ρk ) through (35).

15 16

φ (ρkd )−φ (ρk ) φ (ρk ) Calculate the sensitivity of ρ d : = . ρd δd If (φ (ρkd ) < φ ∗ ), update the optimal solution set (ρ∗ , Y∗ , f∗ , φ ∗ ).} ∗ Calculate step size πk : πk = λk φ (ρk )−ς φ2

17 18

Calculate SR ρ˜k+1 = ρk − πk ∇φ (ρk ) Project ρ˜k+1 on to the feasible space of ρ (0-1), and obtain the new SR ρk+1 . Go to Step 2.}

14

∇φ (ρk )

The step size π k is determined in the same way as in Lo et al (2013) where λk , ς are parameters. From Steps 9 to 15, we change the elements in the vector ρk one at a time until the line frequency fk changes while maintaining the same RTL alignment Yk . It is analogous to the neighborhood searching method around ρk to find a better solution as reflected in Step 15. In summary, the proposed two-phase method improves the existing gradient method by utilizing the dual solutions, and integrates the neighborhood searching method to speed up the solution algorithm. Due to the inherent non-convexity of this problem, the final solution obtained depends on the choice of the starting point ρ, and that global optimality of the solution cannot be guaranteed. The optimal SR obtained from the two phase model is, therefore, a local optimal solution. The performance of the two-phase method will be demonstrated in the case study. 2.4. Stochastic formulation with user equilibrium flows 2.4.1. Model formulation This section formulates the TNDP with stochastic demand under user equilibrium (UE) constraints. From the model formulation standpoint, it can be represented by the combination of stochastic program and bi-level optimization (Patriksson, 2008), which is referred to as stochastic mathematical program with equilibrium constraints (SMPEC). However, the hard capacity constraints of RTL make it inappropriate to simply represent the UE conditions by equilibrium constraints. This issue brings up the need to reformulate the SMPEC with hard capacity constraints. The notion of SR, which separates the decisions over RTL design and passenger assignment into two phases, allows the inclusion of UE consideration or passenger mode-route choice preferences in a convenient manner. Following the framework for SO flows, we separate the formulation into two phases. Phase-1 optimizes the alignment and frequencies of RTL to cover the demand up to ρ. Phase-2 determines the flexible services deployment and UE passenger flows with RTL services determined in phase-1. UE passenger flows under hard capacity constraints and stochastic demand are incorporated in phase-2. As the phase-1 problem here is the same as P1 depicted earlier, in the following we only present the modified phase-2 problem and its solution algorithm catering for UE flows. With flexible services provided along the shortest paths, passengers on the congested RTL services will be diverted to flexible services until the entire system reaches user equilibrium. The modal split between RTL and flexible services is simplified by treating flexible service as the preferred choice when available and hence ρ defines the modal split ratio. We can extend the study to include UE flows across the two modes with preference considerations in the future, but this is out of the scope of this paper. Due to the intrinsic day-to-day demand variations, it is not sensible to find the long-term UE flow pattern. Instead, in this paper, we focus on studying the short-term UE on rapid transit services such that on each day, no passenger taking RTL can decrease his travel time by unilaterally changing routes. Similar to P0-UE, overflow delay for transit line segments operating at capacity is introduced. Similar to the SO formulation, the phase-2 problem endeavors to minimize the expected flexible service and passenger cost through determining the optimal flexible services deployment for each particular demand realization, as expressed below:

Q (ρ, Y, f ) =

 e∈E

pe Qe =



pe (θe + te∗ (θe ) + g∗e (θe ))

(37)

e∈E

where Qe is the phase-2 cost for scenario e, θ e is the total flexible services cost, te∗ (θe ) and g∗e (θe ) are the passenger travel time cost (including in-vehicle travel time, waiting time and transfer penalty) and overflow delay cost, which are functions of θ e . The way to obtain these three cost components is delineated as follows. Given the total flexible services cost fixed at θ e , a linear program problem is formulated to obtain the UE passenger flows in the capacitated RTL through minimizing the passenger travel time cost. The overflow delay of RTL line segment  d,e r r (ir , jr ) is demarcated by the negative Lagrange multiplier −βir,e of the link capacity constraint d∈D Xi j ≤ C (Yi j + Y ji ) f r j r r

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

167

(An and Lo, 2014b). In Section 2.4.2, we will prove that optimizing P3 will produce the UE passenger flows. We then search for the optimal flexible services provision θe∗ so as to minimize the combined phase-2 cost. (P3) Phase-2:

min te (θe ) =



d,e 5 ci4j Xin− vehicle +c

r∈R d∈D

i j∈A d∈D

Subject to:

  1  1 d,e d,e τ + ∗ Xtrans Xboard + c5 f er ∗ 2 fr 2 fr

Xid,e − ir

u∈R∪{0}







Xid,e = ri

u∈R∪{0}

(38)

⎧ d ⎨Be − Zed ,

if i is the origin of OD d d d Z − B , if i is the destination of OD d , ∀i ∈ N, d ∈ D e 0 0 ⎩ e r∈R r∈R 0, otherwise     Xid,e − Xid,e + Xid,e − X jd,e = 0, ∀i ∈ N, r ∈ R, d ∈ D r iu u ir r jr r ir



∀e ∈ E

r∈R d∈D

j∈N



(40)

j∈N

∀i j ∈ A, r ∈ R

Xid,e ≤ C Yirj + Y jir fr , r jr

(39)

(41)

d∈D



cd6 Zed = θe

(42)

d∈D

The objective function te (θ e ) represents the passenger travel time cost consisting of in-vehicle travel time (the first term), waiting time (the second term) and transfer time (the third term). Intuitively, flexible services provide alternatives to overcrowded transit lines and eventually reduce the demand for RTL. After solving P3, the dual solution −βir,e can be obtained j simultaneously. The total overflow delay cost is:

ge (θe ) =





−β

i j∈A r∈R

r,e ij



Xid,e r jr

(43)

d

For a specific demand realization scenario e, by fixing the flexible services provision at the level θ e , we can calculate the passenger travel time cost te∗ (θe ) and the passenger overflow delay cost g∗e (θe ) under UE through P3. The superscript ∗ represents the optimal value of P3. As will be shown in the solution procedure, we develop a line search procedure over θ e to minimize the expected phase-2 cost as expressed in (37). Combining the two phases, the overall objective is to minimize the system cost:

min φ (ρ ) = J (ρ ) + Q¯ (ρ, Y, f ) ρ,θ

(44)

Given a SR ρ, the first term in (44) represents the RTL service cost, which is obtained by solving P1. Upon determining the RTL services, the second term, which stands for the expected phase-2 cost, is solved in P3 and (43) with the flexible services cost fixed at θ = {θe }. We aim to determine the optimal SR ρ∗ and flexible services provision θ ∗ such that the overall system cost φ (ρ) is minimized. We modify the proposed two-phase solution procedure for the SO problem for the inclusion of UE flows, as will be discussed in Section 2.4.3. 2.4.2. Equivalence between P3 and UE assignment By employing the Kuhn–Tucker conditions, following the procedure as described in Bell (1995), Meng et al. (2008) and An and Lo (2014b), we prove that minimizing P3 generates the UE passenger flows under hard capacity constraints, such that passengers on the same OD have the same travel cost. We should note that the inclusion of flexible services Z in P3 distinguishes it from the model described in Bell (1995) and Meng et al. (2008). In addition, the flexible services Z do not bear the same delay as transit services, which is different from the model discussed in An and Lo (2014b). Proposition 2. P3 yields the UE flow pattern under capacity constraint, with the negative Lagrange multiplier associated with the link capacity constraint representing the corresponding passenger overflow delay. Proof. See Appendix.  2.4.3. Solution procedure For the case with SO passenger flows, we only need to optimize ρ, which acts as a given parameter in phase-1. For the case with UE flows, besides ρ in phase-1, we have another variable to be optimized, namely, the cost of flexible services θ e , ∀e, in phase-2. To develop an efficient line search over the two variables ρ and θ e , ∀e, the lower and upper bounds that confine the optimal solution are calculated first. It is trivial that [0, 1] is the feasible region for ρ. As for the flexible services costs θ e in each demand scenario e, we utilize one LP P4 to find its lower bound θ e, min and construct the condition under which the upper bound θ e, max is reached. The bounds form the basis of the line search procedure for the optimal flexible services cost θ e, opt so as to minimize the combined flexible services and passenger cost.

168

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

Note that flexible services are activated only for demand overflow that exceeds the capacity of RTL. The lowest bound cost is to provide just enough flexible services to serve all the realized demands without regard to passenger time cost. In P4, for each demand scenario e, we minimize the flexible services cost with the RTL alignment (Y, f) fixed in P1: (P4)

min θe,min = he ,Ze



cd6 Zed

d∈D

Subject to: (39)–(41) The minimum required flexible services θ e, min determined in P4 to barely fulfill the demand constraints may not be optimal, since it may cause a very high delay cost as many transit line segments may be operating at capacity. Any θ e < θ e, min , however, will make P3 infeasible since certain demand constraints will be violated. Thus, we have θ e, opt ≥ θ e, min . On the other hand, the most expensive way to supply flexible services is such that all passengers are able to take the shortest paths in the RTL network, which will be achieved when the overflow delay is zero, namely, βirj = 0,∀i j ∈ A, r ∈ R. This relationship is proved by contradiction in Lemma 1 as a preparation for Proposition 3. Lemma 1. When the overflow delays are zero for all RTL line segments, βirj = 0,∀i j ∈ A, r ∈ R, all RTL passengers must take the shortest path corresponding to each OD pair. Proof. See Appendix.  Proposition 3. When the overflow delays for all RTL line segments are zero, βirj = 0,∀i j ∈ A, r ∈ R, the required flexible services cost θ e, max is the upper bound of θ e, opt . Proof. See Appendix.  After obtaining the lower bound θ e, min by P4 and the stopping criteria when θ e, max is reached as stated in Proposition 3, we are ready to search for θ e, opt , so that the phase-2 cost is minimized. For a given θ e ∈ [θ e, min , θ e, max ], and with the RTL deployment fixed by P1, we can obtain the passenger cost te∗ (θe ) + g∗e (θe ) by solving P3 and (43). The solution algorithm for the UE formulation can take advantage of the solution procedure developed for the SO case by adding a searching procedure for the optimal flexible services cost θ e, opt . Given an initial ρ, we determine the optimal {Y, f} and the corresponding RTL services cost through solving P1. Then this optimal {Y, f} is utilized as input parameters to P3 and P4. Subsequently, demand scenarios are generated by Latin Hypercube Sampling method to simulate the day to day demand variation and assess the additional flexible service cost. For each scenario, we first calculate the minimum flexible services cost θ e, min by P3. In the context of the ferry network design proposed by An and Lo (2014b), the flexible services provision level is increased by ω each time starting from θ e, min until θ e, max is reached. A smaller step size ω yields a more accurate solution but requires solving the phase-2 problem more times, which is tedious as the phase-2 problem involves a number of LPs. To deal with this challenge, in this study, we develop a two-phase solution algorithm by taking advantage of the special zig-zag shape of the phase-2 cost function. As will be proved in Proposition 4, the phase-2 cost takes on a zig-zag shape and is a non-convex function with respect to θ e . The phase-2 cost function Q(θ e ) is constituted of three cost components, flexible services cost θ e , passenger travel time cost te∗ (θe ) and the overflow delay cost g∗e (θe ). Flexible services cost θ e is obviously an increasing function of θ e itself. Lemma 2 and Lemma 3 analyze the properties of cost components te∗ (θe ) and g∗e (θe ). Lemma 2. The optimal objective function value of P3, te∗ (θe ), is a continuous piecewise linear convex decreasing function of the right hand side value θ e . Proof. See Appendix.  Lemma 3. g∗e (θe ) is a discontinuous, piecewise linear non-convex function with respect to θ e , composed of piecewise constant segments, with each constant segment corresponding to a feasibility range of θ e in which the optimal basis of P3 remains unchanged, and breakpoints, with each breakpoint corresponding to the critical point of θ e that the optimal basis of P3 changes. Proof. See Appendix.  Proposition 4. . The phase-2 cost Q(θ e ) is a discontinuous piecewise increasing function with respect to θ e , composing of piecewise increasing segments and breakpoints, segmented in the same way as g∗e (θe ). The minimum Qe (θ e ) is achieved at one of the drop points, which can be identified by the feasibility ranges. Proof. See Appendix.  In contrast to the condition that the overflow delay g∗e (θe ) is an increasing linear function in the feasibility ranges in An and Lo (2014b), it is a constant in this paper as the overflow delay of flexible services is excluded. Combining the three cost components, the phase-2 cost function is also composed of piecewise linear segments and break points. Each increasing linear segment corresponds to a feasibility range of θ e in which the optimal basis of P3 remains unchanged. Each breakpoint corresponds to a sudden drop in the phase-2 cost where the optimal basis of P3 changes. The number of feasibility ranges depends on the number of vertices in the dual feasible region of P3. This result has important implications on constructing

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

169

3000 2500

Cost

2000 1500

phase-2 cost overflow delay cost flexible services cost passenger travel time cost

1000 500 0 1000

1500

2000

2500

Flexible services cost Fig. 2. The shapes of the cost components as a function of θ e . (For interpretation of the references to color in this figure, the reader is referred to the web version of this article).

the solution algorithm, as the optimal solution must be achieved at one of the drop points where the optimal basis changes. The feasibility range of θ e can be readily obtained after solving P3 and typically provided by many LP solvers, such as CPLEX. It enables us to jump to the drop points directly as they are also the endpoints of the feasibility ranges, which avoids enumerating all points in [θ e, min , θ e, max ]. The one with the lowest phase-2 cost is chosen as the optimal solution. This method not only relieves the otherwise cumbersome enumeration process in the region of [θ e, min , θ e, max ], but also enables us to obtain the optimal solution instead of a near optimal solution subject to the searching step size ω. The number of iterations required depends on the number of vertexes in the dual feasible region of P3. This critical observation greatly improves the algorithm efficiency in searching for the optimal θ e, opt in the phase-2 problem. Repeatedly solving P3, P4 and the searching process within [θ e, min , θ e, max ] for each demand scenario e, we obtain the expected flexible services cost and expected passenger cost for a given {Y, f} under the short-term UE passenger flows. In the following, we delineate the modified phase-2 solution procedure. For a given {Y, f}, repeat the following procedure for each demand realization e = 1 . . . m Step 1. Solve P4 to get the minimum flexible cost θ e, min and set θe = θe,min . Step 2. For fixed θ e , solve P3 to get passenger cost te∗ (θe ), phase-2 cost Qe = θe + te∗ (θe ) + g∗e (θe ), and the feasibility range [θe , θe ] of θ e . Step 3. If Qe < Qe∗ then set Qe∗ = Qe . Step 4. If g∗e (θe ) = 0, stop and take the optimal phase-2 cost as Qe∗ . Otherwise set θe = θe + ε , go to step 2.  The expected phase-2 cost is calculated as the sample average of demand scenarios e ∈ E, Q (ρ, Y ) = e∈E pe Qe∗ . Note that ɛ in Step 4 is a small positive number to push θ e into the next adjacent feasibility range. The searching method for θ e, opt is justified due to the non-convexity of the phase-2 cost Qe with respect to flexible services provision θ e . To illustrate the shapes of the phase-2 cost components functions with respect to θ e , we plot out the cost components for one demand scenario arbitrarily selected from our case study, as shown in Fig. 2. The x-axis stands for the amount of flexible services provided θ e . The black dotted line represents the flexible cost θ e with a slope of 1. The passenger travel time cost te∗ (θe ), represented by the green dotted dashed line, is a linearly decreasing function. The overflow delay g∗e (θe ), represented by the red dashed line, is a non-convex function composed of flat segments and break points where sudden drop occurs. The flat segment represents the feasibility ranges [θe , θe ] of θ e . Combining the three cost component, we plot the phase-2 cost in the blue solid line. It is a zig-zag and non-convex function with piecewise increasing linear segments. As can be observed from the shape of Qe (θ e ), the optimal level of flexible services provision θe∗ can only be achieved at the breakpoints, where the optimal basis of P3 changes. The solution algorithm makes use of this observation and only calculates the phase-2 cost at the break points through detecting the boundaries θe , θe of the feasible ranges for each specific θ e . 3. Case study 3.1. An illustrative network The same network and settings as in An and Lo (2014a) are adopted as shown in Fig. 3. It contains 10 nodes and 38 directed links with their distances lij shown by the side of the links. Up to three rapid transit lines are constructed, i.e. Rmax = 3, with the operating frequency bounded between [ fmin , fmax ] = [3, 20] transit units (TU) per hour. The capacity of a transit unit is C = 80 passengers/TU. The average hourly demands bd of the 9 OD pairs are presented in Table 2. The random demands follow the lognormal distribution Bd ∼ log N(bd , σ 2 ) with the expectation of bd and variance of σ 2 = 4∗ bd . The

170

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

9

5

10

8

9

4

7 9

7

3

3

9

10

2

5

12

8

10

1

6

6

3

13

4

5 3

6 16

7 3

Link cost Fig. 3. Network setting.

unit operating cost per vehicle per unit length is c1 = 1 and the unit length construction cost is c2 = 10. The unit station construction cost is c3 = 1. The unit cost per passenger per length for flexible services is set at four times the RTL operating cost, i.e. c6 = 4/C. The on-vehicle passenger value of time is c4 = 0.01. The waiting cost per unit time for passengers is c5 = 0.02 and the transfer penalty is τ = 1. The sample size adopted in phase-2 is m = 100. The initial value of ρ is 0.5 and the two-phase algorithm starts the searching from the deterministic solution. The proposed two-phase method was programmed in C++ and the MILP was solved by the IBM CPLEX software Academic version 12.3 and run on Intel i7-3370 computer with 16GB RAM, 64 bit. In this paper, frequencies are assumed to be continuous for computational efficiency. However, in practice, frequencies in real numbers may not be practical. To mitigate this problem, an enumeration approach is developed. We assume each RTL headway is an integer multiple of 30 s, 60×2 and and each optimal frequency obtained from the robust or stochastic model has two possible candidates: fr∗ = 60×2/ f ∗ r

60×2 fr∗ = 60×2/ , where ·  and  ·  are the floor and ceiling functions, respectively. If a total of Rmax RTL routes are to fr∗  be constructed, then there will be up to 2Rmax frequency combinations to be evaluated to get the system cost. The solution with the lowest system cost is used to represent the optimal solution for real world applications. Table 1 shows the optimal results under SO and UE passenger flows by the proposed two-phase method. To evaluate the performances of the robust model and two-stage stochastic models, we determine the optimal solutions of the traditional robust approach with stochastic demand for three uncertainty sets, ranging from 1 × σ , 2 × σ and 3 × σ under SO passenger flows, and present the results in Table 1. Given the RTL alignments and frequencies obtained by solving the robust formulation P0-SO , the expected passenger and flexible service costs of the robust solution are determined by solving the phase-2 problem, namely P2 for the SO case, and P3 and P4 for the UE case, which simulates or evaluates the expected performances of the RTL solutions of the robust formulation under SO and UE, respectively. Hence the cost components shown in Tables 1 and 3 reflect the actual costs and thus are comparable among the traditional robust model and the two-phase model. Statistically, 3 × σ can cover up to 95% of demand realizations for most stochastic distributions. The alignments and frequencies for the three lines, cost components as well as computation times are listed in each column. The service reliability is the arithmetic average of SR over all OD pairs, namely, the arithmetic average of elements in ρ. The designed SR, as shown in the third column of Table 1, reflects the optimal ρ as determined with the two-phase method P1, and the input based on the uncertainty bound of the robust model P0-SO . It indicates that the RTL network is able to handle the worst scenario when all the OD demand realizations simultaneously reach their maximum values, as limited by ρ. However, in reality not all the demands would reach their maximum values at the same time in RTL operations. Hence the actual demand that can be carried by RTL is higher than the designed value. The actual SR as shown in the fourth column of Table 1 reveals the average RTL patronage. It is calculated as the average ratio of passengers taking RTL to the expected demand over all OD pairs. We observe that the solutions obtained by the proposed two-phase method possess the lowest total costs in their corresponding categories: UE and SO. Under SO passenger flows, the two-phase method deploys the least RTL services and relies much on flexible services. Under SO, as the breadth of the uncertainty set increases from 0 × σ to 3 × σ , the RTL cost increases from 1156 to 1456 to cater for the higher SR; the flexible services cost almost decreases to zero under the 3 × σ scenario with a SR of up to 0.998. The two-phase method and the cases of 0 × σ ,1 × σ , 2 × σ generate different alignments and frequencies, many with track sharing links. The track sharing lines allow effective use of RTL tracks, hence reduce the line and station construction cost. The optimal SR obtained by the two-phase method is [0.73, 0.87, 0.66, 0.54, 0.64, 0.96, 0.64, 0.83, 0.70] with an average of 0.72. For this case study, the two-phase method (robust formulation) requires 30 (2– 5) and 51 (8–12) s of computational times for SO and UE flows, respectively.

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

171

Table 1 The line alignments and cost components for solutions under SO and UE passenger flows.

SO

UE

a b c

Service reliability

Line

Freq.

Company cost

Expected passenger cost

Expected

Comp.

Designed

alignment

(veh/h)

RTL

RTL

Actual

2-phase method

0.72

0.960

0 × σb

0.50

0.948

1 × σb

0.84

0.992

2 × σb

0.97

0.999

3 × σb

0.998

1.00

2-phase method

0.79

0.978

0 × σb

0.50

0.935

1 × σb

0.84

0.879

2 × σb

0.97

0.979

3 × σb

0.998

0.998

4–5–6–1–2–8–7 5–6–1–2–8–9 Z3–2–10 5–6–1–2–8–9 3–1–2–10 4–1–7 8–9–10–3–1–6–5 7–8–9–10–3–4 6–1–2–8 4–5–6–1–2–8–9–10 5–6–1–2–8–7 4–3–10–9–8 9–8–2–1–6–5 4–1–7 3–2–10

10 10 6.7 12 6.3 10 7.5 11 11 6.7 13.3 13.3 17 15 10

9–8–2–1–6–5 4–1–7 3–2–10 5–6–1–2–8–9 3–1–2–10 4–1–7 8–9–10–3–1–6–5 7–8–9–10–3–4 6–1–2–8 4–5–6–1–2–8–9–10 5–6–1–2–8–7 4–3–10–9–8 9–8–2–1–6–5 4–1–7 3–2–10

13.3 12 7.5 12 6.3 10 7.5 11 11 6.7 13.3 13.3 17 15 10

Flexible

Flexible

total cost

time (s)

a

45

2044

30

1070

180

713|20|16 749

1156

158

677|23|16a 16

40

2070

4

1249

27

833|21|17a 871

6

2153

2

1363

2

836|17|16a 869

0.1

2234

3

1456

1

706|22|11a 739

0

2196

5

1265

56

690|22|14|7c 733

14

2068

51

1156

167

673|23|16|44c 756

42

2121

8

1249

94

797|18|16|74c 905

24

2272

10

1363

14

828|17|16|52c 917

4

2298

10

1456

8

704|22|11|1c 738

2

2204

12

Passenger cost components taking RTL under SO: in-vehicle cost| transfer cost| waiting cost; the bottom one: total passenger cost. n × σ stands for solutions with uncertainty set defined as [bd − n∗ σ , bd + n∗ σ ]. Specifically, n = 0 stands for the deterministic solution. Passenger cost components taking RTL under UE: in-vehicle cost| transfer cost| waiting cost| overflow delay cost; the bottom one: total passenger cost.

After including the overflow delay for the UE assignment, the passenger cost on transit lines increases, making passenger cost larger. It is sensible that the company decides to decrease the overflow delay by increasing the frequencies of the transit lines or calling upon more flexible services to handle demand uncertainty. This explains the big increase in transit operating costs as compared with the SO counterpart solution at the optimal SR. Since the transit lines can handle more passengers with higher frequencies, the need for flexible services declines. It is the same case for passenger cost of flexible services. The optimal SR is [0.86, 0.86, 0.82, 0.68, 0.59, 0.93, 0.86, 0.70, 0.87] with an average of 0.79. It is a little higher than that of the SO solution, indicating that the SR is enhanced in this case if passenger overflow delay is taken into account. Under UE, the two-phase method generates a different RTL alignment to bring down the overflow delay to 7. For the traditional robust method with a predefined level of robustness, the overflow delay is up to 74. The cost difference between stochastic solutions obtained by the two-phase method and robust solutions is up to 190 under SO and 230 under UE. It shows that the system can save up to 9% under SO and can save as high as 11% under UE, should the demand distribution information be made available. This value of information could be higher for a more dispersed distribution with a larger variance. The solutions obtained by the two-phase method are selected to represent the SO and UE solutions. The alignments of the three rapid transit lines are plotted in Fig. 4. We can observe that the resultant three RTLs cover all ten nodes and are interconnected by transfer stations (labeled in grey) into a fully connected network under both UE and SO conditions. However, the RTL network configuration under SO and UE are starkly different. In particular, the RTL network under SO has a total transit length of 48 with 9 links and five transfer stations, nodes 1, 2, 5, 6 and 8. The RTL network under UE consists of a total length of 60 with 9 links and two nodes 1 and 2 are set for transfer stations. This result indicates that including the UE conditions compels the system to construct the RTL network with less links but longer distances, such as links 1–7 and 1–4. Although long links are more expensive, and thus increase the company’s cost, it enhances the directness of the RTL, which is beneficial to passengers. Specifically, it increases the company’s construction and operating cost from 1070 to 1265 but substantially reduces the passenger cost and flexible services cost from 974 (= 180 + 749 + 45) to 803 (= 56 + 733 + 14). The resultant average traffic flow patterns under SO and UE passenger flows are shown in Table 2a and b. The RTL services carry most of the passengers, around 94% of all travel demand for the SO solution and 98% for the UE solution with flexible services serving as a supplementary role as expected. The last column of the two tables shows the distance of shortest path given by flexible services, and the third column from right gives the path distance on the RTL. The subscripts

172

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

Fig. 4. Resultant RTL configurations and passenger flow.

Table 2 OD demand and resultant passenger flow. OD pair

OD demand

Path by RTL

Expected RTL patronage

Transfer

Distance

(a) Solutions under SO 2→10 200 23 –103 3→2 150 33 –23 4→7 800 41 –51 –61 –11 –21 –81 –71 5→8 350 51 –61 –11 –21 –81 52 –62 –12 –22 –82 6→9 600 62 –12 –22 –82 –92 7→6 250 71 –81 –21 –11 –61 8→3 400 82 –22 –23 –33 9→4 450 92 –82 –81 –21 –11 –61 –51 –41 92 –82 –22 –12 –62 –52 –51 –41 10→5 500 103 –23 –22 –12 –62 –52

200 150 743 32 318 412 250 395 429 21 474

0 0 0 0 0 0 1 1 1 1 1

8 10 27 17 17 17 18 16 26 26 19

(b) Solutions under UE 2→10 200 23 –103 3→2 150 33 –23 4→7 800 42 –12 –72 5→8 350 51 –61 –11 –21 –81 6→9 600 61 –11 –21 –81 –91 7→6 250 72 –12 –11 –61 8→3 400 81 –21 –23 –33 9→4 450 91 –81 –21 –11 –12 –42 10→5 500 103 –23 –21 –11 –61 –51

200 150 792 350 567 250 399 422 392

0 0 0 0 0 1 1 1 1

8 10 22 17 17 16 16 23 19

Expected flexible patronage

Distance

0 0 57 0

8 10 22 17

188 0 5 0

15 16 16 21

26

19

0 0 8 0 33 0 1 28 8

8 10 22 17 15 16 16 21 19

A2 –B2 –B3 –C3 stands for a path from station A to C via station B, with station B serving as transfer station from line 2 to line 3.

in the column of path by RTL represent line index. In the RTL network under SO, passengers on one OD pair, such as OD pairs 4–7, 9–4 etc., may be assigned to different paths with significant longer distance as compared with that on flexible services. Although 5 stations can serve as transfer stations between line 1 and line 2, most transfers happen at the diverging stations of the two lines, i.e. stations 5 and 8. For the UE solution, in contrast, there are no alternative paths for all OD pairs and the path lengths are much shorter, with most of them having the same shortest distances as those of flexible services. Specifically, for the current RTL design under UE, the used path is the unique feasible path for each OD pair. Constructing more lines together with higher frequencies greatly decreases the chance of having the transit lines operating at capacity, wherein the overflow delay occurs. It thus prevents passengers from long travel distance on the RTL, which explains why the overflow delay for the final design is only 7. Under both SO and UE, most passengers can take direct services and the number of transfer is at most 1. Multiple transfers do not happen under SO, although allowed, possibly due to the high transfer penalty. The UE consideration does substantially influence the RTL design and the total passenger cost.

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

173

Fig. 5. Resultant RTL configurations and frequencies of the robust formulation.

Table 3 The cost components, level of robustness and computation time for solutions under SO and UE passenger flows. Level of robustness

Company cost

Expected passenger cost

Expected

Comp.

Designed

Actual

RTL

Flexible

RTL

Flexible

total cost

time (s)

SO

Two-phase method 0×σ 1×σ 2×σ

0.54 0.50 0.84 0.97

0.88 0.91 0.99 1.00

1765 1936 2384 2965

446 275 0.01 0

1313 1413 1437 1388

112 69 0.01 0

3636 3693 3821 4353

211 7200 7200 7200

UE

Two-phase method 0×σ 1×σ 2×σ

0.55 0.50 0.84 0.97

0.91 0.87 0.99 0.99

1892 1936 2384 2965

324 438 28 18

1517 1468 1484 1392

81 110 7 5

3814 3952 3903 4380

5005 7200 7200 7200

3.2. Sioux falls network To demonstrate the applicability of the proposed model for a larger network, we apply it to the Sioux-Falls network with 24 nodes and 76 directed links (Bar-Gera, 2009) as shown in Fig. 5. A total of 226 OD pairs with demands larger than 500 passengers/h are selected to represent the passenger demand in the network. The capacity of a transit unit (TU) is set to be C = 1600 passengers/h and the maximum frequency fmax is set at 20 TU/h, offering a maximum line capacity of 32,000 passengers/h. The in-vehicle value of time is c4 = 0.0005. The unit passenger waiting cost is c5 = 0.001 and the transfer penalty is τ = 1. The other parameters follow the settings in the previous example unless otherwise specified. For this larger case study, we set the computation time limit to be 2 h. The traditional robust method is not able to identify the optimal solution within this time limit. A residue gap of 3% cannot be cleared even after 2 h of computation. Actually, we run the 0 × σ case for 5 h, but a residual gap of 1.7% still cannot be cleared, and the best feasible solution does not improve in the last three hours of computation. In Table 3, the computation time for the robust formulation in the right-most column only reports the 2-h time limit, not the actual time of achieving optimality. In the proposed two-phase method, we modify the phase-1 problem by adding the flexible service cost, as shown in P1 . Since flexible service deployment Zd is a real variable, it somehow relaxes the strict equality constraint (46), which is able to reduce the computation time of phase-1 substantially from more than 2 h to around 30 s. Under SO, the computation time of the two-phase method is 211 s, significantly less than 2 h. The computation time for the UE case increases to 5005 s due to two reasons: 1) the phase-2 problem under UE is calculated repeatedly to find the optimal flexible service cost θ e ; 2) the gradient calculation over ρ must go through the tedious solution process of phase-2 repetitively, which is evaded for the SO solution. In future studies, it would be beneficial to approximate the gradient of ρ under UE to improve the algorithm efficiency. Comparing the computational performances for the small illustrative network and the larger Sioux Falls network, it appears that the advantage of the two-phase method is more pronounced for the larger network.

174

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

Fig. 6. Resultant RTL configurations and passenger flows for the two-phase method.

(P1 ) Modified Phase-1:

min f,Y



ci1j yrij +

r∈R i j∈A



ci2jY i j +

i j∈A



c3Wi +

i∈N



d ci4j Xin− vehicle +

i j∈A d∈D



cd6 Z d

(45)

d∈D

subject to (2)–(18), (20)–(21) and (26)

 r∈R

Xid0 ir −

 r∈R

Xidr i0 =

⎧ d ⎨B¯ − Z d , Z d − B¯ d , ⎩ 0,

if i is the origin of OD d if i is the destination of OD d , otherwise

∀i ∈ N, d ∈ D

(46)

The line alignments and frequencies obtained by the two-phase method are shown in Fig. 6, while the solutions obtained by the traditional robust optimization with different breadths n are shown in Fig. 5. Fig. 5 shows that all nodes except node 3 are constructed as RTL stations since no demand is generated or ended at node 3. In the solutions obtained by the twophase model as shown in Fig. 6, nodes 1, 2, 3 and 6 are not coved under SO, and nodes 1, 2 and 3 are not coved in UE. Due to low demand at these nodes, it is not economical to stretch the RTL to these far-end nodes. Instead, flexible services are employed to carry the demands not covered by RTL. Traditional robust method fails to distinguish the nodes with low demand and treats all nodes equally. The solution includes multiple lines sharing the same rail track, such as links 10–1, 15–22, etc. in Fig. 5a. Since symmetric OD demands are adopted in the Sioux Falls network, the passenger flows on the two directions of a link are similar with less than 0.1% difference. The average RTL link flows of the two directions are calculated and presented along each link in Fig. 6. Note that the values are scaled down by 100 for better presentation. The RTL alignments are substantially different under SO and UE as determined by the two-phase method. In particular, 5 track-sharing links are observed under SO, whereas only 2 track-sharing links are observed under UE. Based on the experiments, cyclic route solutions never occur in the small illustrative example, while cyclic route solutions do occur in the Sioux Falls network. To rectify this problem, after identifying cyclic routes through the SO algorithm, we add the cycle breaking constraint (18) for the observed cycles in the UE and traditional robust model. The SO model can identify cyclic solutions in about 30 s of computational time for the Sioux Falls network. The added cycle breaking constraints, about 10 for this network, can evade the cyclic route problem, and hence reduces the computation times for solutions. 3.3. Applications to other network configurations To further demonstrate the applicability of the two-phase models, this section conducts a series of sensitivity analyzes on four network configurations as shown in Fig. 7. Link distance is shown along each link. We vary the maximum route number allowed Rmax from 1 to 3 as shown in Table 4, and high and low demand scenarios are generated. Matrix b represents the expected OD demands of the grid network while b is used in the other three networks. In the high demand scenarios, the expected demand is doubled. The other parameters follow the settings in the Sioux Falls Network.

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

175

Fig. 7. Network settings. (a) Grid network (Guan et al., 2006). (b) Triangle network (Marin and Jaramillo, 2008). (c) Bone network. (d) Star network. Table 4 Computation results for the four networks. Network

Grid Grid Grid Grid Grid Grid Triangle Triangle Triangle Triangle Triangle Triangle Bone Bone Bone Bone Bone Bone Star Star Star Star Star Star

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ b = 100∗ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎜ ⎜ ⎜ ⎜ ∗⎜ b = 100 ⎜ ⎜ ⎜ ⎜ ⎝

Demand

Low Low Low High High High Low Low Low High High High Low Low Low High High High Low Low Low High High High

Rmax

1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

SO

UE

Same line

RTL cost

Total cost

LOR

Comp. Time (s)

RTL cost

Total cost

LOR

Comp. time (s)

align. under SO&UE

887 1134 1072 1629 1693 1695 1358 1133 1080 1536 1535 1555 520 1128 1257 698 1353 1572 407 839 1358 580 1230 1681

2181 1902 1822 4224 3127 3031 2647 1998 1916 4283 3390 3274 4002 2932 2277 7474 5443 4841 4060 3208 2374 7646 5919 4611

0.81 0.95 0.97 0.70 0.96 0.98 0.89 0.96 0.97 0.78 0.93 0.97 0.20 0.71 0.96 0.20 0.57 0.73 0.22 0.51 0.97 0.21 0.48 0.81

1 4 25 1 7 105 6 16 35 3 30 375 2 2 3 2 3 3 1 3 3 1 1 2

896 1192 1142 1754 1524 1710 1430 1205 1110 1582 1725 1599 505 1147 1279 710 1444 1592 418 932 1400 572 1180 1673

2192 1931 1843 4279 3300 3109 2714 2051 1961 4310 3557 3406 4021 2964 2309 7499 5586 4990 4108 3283 2427 7706 5991 4766

0.78 0.94 0.97 0.64 0.93 0.99 0.93 0.95 0.96 0.73 0.97 0.98 0.16 0.71 0.96 0.19 0.63 0.74 0.20 0.50 0.96 0.19 0.47 0.78

3 5 30 10 30 330 131 144 299 2 651 599 4 61 535 4 45 54 1 43 1 5 15 290

Y Y N Y Y N N N Y Y N Y Y Y Y Y N N Y Y Y Y Y N

− 0 30 30 30 0 0 0 0

0 − 22 0 22 0 0 0 0

30 22 − 30 52 0 0 0 0

30 0 30 − 40 88 0 0 0

30 22 52 40 − 0 52 0 52

0 0 0 88 0 − 42 42 0

0 0 0 0 52 42 − 42 52

0 0 0 0 0 42 42 − 0

0 0 0 0 52 0 52 0 −

− 9 26 19 13 12 13 8 11

9 − 14 26 7 18 3 6 12

26 14 − 30 24 8 15 12 5

19 26 30 − 22 16 25 21 23

13 7 24 22 − 20 16 22 21

12 18 8 16 20 − 16 14 12

13 3 15 25 16 16 − 11 11

8 6 12 21 22 14 11 − 4

11 12 5 23 21 12 11 4 −

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

176

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

We can observe that the total system cost always drops when Rmax increases under both low and high random demands and SO/UE conditions. With more RTL lines allowed to be constructed, the RTL lines can align in a more economical way such that the passenger patronage to each segment of the line does not vary too much. A balanced RTL utilization leads to a lower total system cost, higher RTL level of robustness (LOR) and sometimes even lower RTL construction cost. When demand is doubled, longer lines with higher frequencies are constructed to handle the demand increase. However, the RTL level of robustness could drop at high demand level probably because some lines have reached the maximum frequency limit and are no longer able to take extra passengers. The UE passenger flow condition can change the line alignment as indicated in the last column of Table 4. Specifically, the triangle network is more prone to be affected by UE while the star network is not, since the triangle network has more path alternatives for each OD pair than those of the star network. As far as computational effort is concerned, all cases can be solved within 12 minutes, demonstrating the efficiency of the two-phase solution approach. 4. Conclusion In this paper, we proposed a two-phase stochastic program and a service reliability based solution algorithm for rapid transit network design under demand uncertainty. Flexible services were introduced to carry the demand overflow, which could be interpreted as the expected cost incurred by outsourcing to a third party to carry the demand overflow, or simply the cost of demand loss due to inability to carry the demand overflow by RTL. By minimizing the combined expected cost of the regular and flexible services, the local optimal service reliability was obtained by the proposed two-phase solution algorithm. The two-phase solution procedure developed integrates an improved gradient algorithm with a neighborhood search, which substantially reduces the computation time as compared with the gradient method proposed by An and Lo (2014a). In applying the two-phase method for the Sioux Falls network, the results were promising. In addition, we compared the system performances between the two-phase method and the robust approach with different demand uncertainty sets. The results showed that substantial costs could be saved if the service reliability could be optimally determined rather than defined a priori. The difference in cost between the two approaches provides an estimate of the value of gathering further information on the demand distribution. Lastly, the results showed that the inclusion of UE flows markedly changed the RTL alignments and cost components. In the current study, we assumed flexible services as an alternative mode, with the advantage of shorter-distant or more direct routes to passengers. In practice, over-deployment of flexible services may cause congestion on the roadway network, hence may substantially increase the overall network travel time. One possible future extension is to take into account the roadway congestion caused by flexible services, which will refine the RTL design. Nevertheless, we believe that the notion of optimal service reliability introduced here will offer a new way to consider transit network design under demand uncertainty. Acknowledgments The study is supported by the Public Policy Research Grant HKUST6002-PPR-11 and General Research Fund #616113 of the Research Grants Council of the HKSAR Government, the Hong Kong Ph.D. Fellowship, and Strategic Research Funding Initiative SRFI11EG15. Appendix Proposition 2. P3 yields a UE flow pattern under capacity constraint, with the negative Lagrange multiplier associated with the link capacity constraint representing the corresponding passenger overflow delay. Proof. In the following proof, the subscript e is omitted for the concerned variables and each variable stands for a particular scenario. Let hdκ be the passenger path flow on path κ for OD pair d and with the path set denoted by Kd . The subscript ij in passenger flow Xidj here is simplified to notate either in-station or out-station links. The link-path conservation equation    can be expressed as d∈D Xidj = d∈D κ ∈K d ηidj,κ hdκ where ηidj,κ is 1 if path κ connecting OD pair d traverses link ij and 0 otherwise. We remove the subscript e for the concerned variables as well. The passenger flow conservation Constraints (39)  and (40) can be rewritten in path flow as Bd − Z d = κ ∈K d hdκ , ∀d ∈ D. The Lagrange for P3 can be written as following:

L

h,Z,α ,β ,γ

=

 i j∈A

ci j

 d∈D

Xidj

+

 d∈D



α

d

d

d

B −Z −

 κ ∈K d

d

hκ

+

 i j∈A



β

r ij

 r

C



Yirj

+ Y jir



fr −

 d



Xidj

+γ

θ−



cd6 Z d

(47)

d∈D

where cij is the general unit traveling cost, equals ci4j if ij is an out-station link, and equals the unit waiting cost

c5 2 fr

or

transfer cost c5 (τ + 21f ) if ij is an in-station link as depicted in (38). α d , βirj , and γ denote the corresponding Lagrange mulr tipliers. Note that βirj is defined for out-station links which have capacity restrictions. The in-station links where passengers

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

177

make transfer or boarding do not have capacity restrictions; hence their corresponding Lagrange multipliers are always zero. The Kuhn–Tucker conditions for P3 can be rewritten as:



hdκ 



ci j ηidj,κ −

i j∈A

ci j ηidj,κ −

i j∈A

βirj ηidj,κ − α d = 0 ∀κ ∈ K d , ∀d ∈ D

(48)

i j∈A

βirj ηidj,κ − α d ≥ 0 ∀κ ∈ K d , ∀d ∈ D

(49)

i j∈A

hdκ ≥ 0



βirj 









∀κ ∈ K d , ∀d ∈ D 



C Yirj + Y jir fr −

r



C Yirj + Y jir fr −

r





(50)

Xidj

=0

∀i j ∈ A, r ∈ R

(51)

d

Xidj ≥ 0

∀i j ∈ A

(52)

d

βirj ≤ 0 ∀i j ∈ A, r ∈ R   Z d −α d − γ cd6 = 0 ∀d ∈ D

(54)

−α d − γ cd6 ≥ 0

(55)

(53)

∀d ∈ D

Z d ≥ 0 ∀d ∈ D  θ− cd6 Z d = 0 d∈D d

Bd − Z =

 κ

hdκ

(56) (57)

∀d ∈ D

(58)

∈K d

Constraints (48)–(58) constitute the complementarity conditions for UE passenger flow: when the passenger path flow hdκ   is greater than 0, then the term ( i j∈A ci j ηidj,κ − i j∈A βirj ηidj,κ ) − α d = 0 or that the path cost for OD pair d is the summation of free flow travel time cost and the overflow delay on all links ij connecting path κ . The over flow delay −βirj is positive if the link operates at capacity, zero otherwise. Inequality (49) ascertains that α d is the lowest cost for OD pair d, as all paths connecting d have a higher or equal cost. Constraints (51)–(53) define the overflow delay term −βirj . According to   d r r (51), whenever βirj is negative, the sum of transit service capacity r C (Yi j + Y ji ) f r equal the passenger link flow d Xi j ,

i.e. the capacity constraint is active. Constraints (54)–(57) state the conditions for flexible services Zd to be deployed. This completes the proof for UE passenger assignment.  Lemma 1. When the overflow delays are zero for all RTL line segments, βirj = 0,∀i j ∈ A, r ∈ R, all RTL passengers must take the shortest path corresponding to each OD pair. Proof. When all RTL line segments have zero overflow delay, this indicates that all line segments have excess capacity or that the capacity constraints for RTL services (41) are not binding. Assume that some passengers from OD d do not take the shortest path κ 0 but a path κ 1 with longer distance. Since extra seats are available on every RTL line segment of path κ 0 , these passenger are able to change routes from κ 1 to κ 0 to further decrease their travel time while not influencing the other passengers. The solution with passenger choosing a longer path κ 1 contradicts with the objective function (38) in which the total passenger travel time cost ought to be minimized; hence it is not the optimal solution of P3 when βirj = 0,∀i j ∈ A, r ∈ R. We therefore conclude that when the overflow delay for every RTL line segment is zero, passengers on RTL services must take the shortest paths.  Proposition 3. When the overflow delays for all RTL line segments are zero, βirj = 0,∀i j ∈ A, r ∈ R, the required flexible services cost θ e, max is the upper bound of θ e, opt . Proof. Let the case be called ξ when βirj = 0,∀i j ∈ A, r ∈ R, and the associated required flexible services cost represented by

θ e, ξ . According to Proposition 2, the Kuhn–Tucker conditions related to path flow hdκ on path κ of OD pair d for minimizing P3 can be rewritten as:



hdκ

 i j∈A

ci j ηidj,κ −



βirj ηidj,κ − α d = 0 ∀κ ∈ K d , ∀d ∈ D.

i j∈A

 If the overflow delays βirj are all zero, the path travel time for OD pair d is α d = i j∈A ci j ηidj,κ , meaning that all passengers on RTL are able to select the shortest paths without delay as proved in Lemma 1. The total passenger cost is then the

178

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

product of the shortest path travel time α d and the amount of passengers on RTL, i.e. te∗ (θe,ξ ) = with condition of zero overflow delays g∗e (θe,ξ ) = 0, the phase-2 cost is written as:

Qe

 d∈D

α d (Bde − Zed ). Together



      6 d   d  d   cd Ze θe,ξ + α ∗ Be − Zed θe,ξ θe,ξ = θe,ξ + te∗ θe,ξ + g∗e θe,ξ = =



α d Bde +

d∈D





d∈D

cd6 − α d ∗ Zed



θe,ξ

d∈D



d∈D

Note that α d is a constant. Zed (θe,ξ ) stands for the flexible services deployment under a total flexible services cost of θ e, ξ . For any θ e > θ e, ξ , passengers will shift from RTL to the flexible services as they provide door to door services with shorter

travel time, leading to Zed (θe ) ≥ Zed (θe,ξ ). Since the unit flexible services cost cd6 (operating and passenger cost combined) is

set to be larger than the passenger travel time cost α d on RTL, we have cd6 − α d ≥ 0. Combining the condition that Zed (θe ) ≥ Zed (θe,ξ ), we can get for any θ e > θ e, ξ that:

Qe (θe ) =



α d Bde +

d∈D





cd6 − α d ∗ Zed (θe ) ≥

d∈D



α d Bde +

d∈D





cd6 − α d ∗ Zed



   θe,ξ = Qe θe,ξ .

d∈D

Or that for any θ e > θ e, ξ , the phase-2 cost follows Qe (θ e ) ≥ Qe (θ e, ξ ). Hence the optimal flexible services cost θ e, opt to minimize the phase-2 cost Qe must be less than θ e, ξ . In other words, the value of θ e, ξ that corresponds to βirj = 0,∀i j ∈ A, r ∈ R is an upper bound, or θe,ξ = θe,max . Combining with the condition that θ e, opt ≥ θ e, min , we have θ e, opt ∈ [θ e, min , θ e, max ].  Lemma 2. The optimal objective value of P3 te∗ (θe ) is a continuous piecewise linear convex decreasing function of the right hand side value θ e . Proof. It is proved that the optimal objective function value of a linear programing problem is a convex function with respect to the right hand side value of the constraint set. More specifically, it is a piecewise linear function (Bertsimas and Tsitsiklis, 1997). Hence t∗(θ e ) is a piecewise linear convex function with respect to the right hand side value θ e of the  6 d 6 d d resource constraint d∈D cd Ze = θe . According to (55), we have γ cd ≤ −α ∀d ∈ D, where α is positive representing the generalized path cost for OD pair d. Thus we have γ cd6 ≤ −α d ≤ 0. With the Lagrange multiplier γ ≤ 0 for the resource  constraint d∈D cd6 Zed = θe , f∗(θ e ) is a decreasing function with respect to θ e . Together with the convexity and piece-wise linearity of f∗(θ e ) , f∗(θ e ) is a piecewise linear convex decreasing function. 

Lemma 3. g∗e (θe ) is a discontinuous, piecewise linear non-convex function with respect to θ e , composed of piecewise constant segments, with each constant segment corresponding to a feasibility range of θ e in which the optimal basis of P3 remains unchanged, and breakpoints, with each breakpoint corresponding to the critical point of θ e that the optimal basis of P3 changes. Proof. The dual feasible region of P3 is a polyhedron with a finite number of vertices which does not change with the resource parameter θ e . Let α, β, γ , whose components defined in (38), are the optimal solution of the dual P3, which must be on the vertices or boundaries of the dual feasible region, and may change from one vertex to the other as θ e varies. We separate the analysis of g∗(θ e ) into two cases: one is along any feasibility range in which the optimal basis of P3 does not change; the other is any breakpoint at which the optimal basis changes. Case 1: First, we consider the range of θ e such that the optimal basis is preserved. i.e. the dual optimal solution α, β, γ does not change with the right hand side θ e . By differentiating the overflow delay g∗(θ e ) with respect to θ e , we obtain:

∂ (g∗ ) ∂ = ∂θe



i j∈A

−βirj

∂θe



d

Xidj



=−

 i j∈A





  ∂ d Xidj r ∂βirj  d X + βi j ∂θe d i j ∂θe

(59) ∂β r

In the feasibility range or when the optimal basis remains, all the Lagrange multipliers are unchanged, that is ∂θi j =0. Ace   cording to constraints (51)–(53) in the proof of Proposition 1, βirj < 0 indicates that r C (Yirj + Y jir ) fr = d Xidj , then we have ∂(



Xidj )

∂



∂(



d d Xi j ) and βirj is always zero. Combining the two conditions, ∂θe ∂θe ∂ (g∗ ) we can obtain ∂θ = 0. The overflow delay g∗(θ e ) is a constant function of θ e in the feasibility range of a particular θ e . The e d

∂θe

=

r r r C (Yi j +Y ji ) f r

= 0. It indicates that the product of

number of feasibility ranges depends on the number of vertices in the dual feasible region of P3. Case 2: now we consider the break points at which the optimal basis of P3 changes.

According to Proposition 3, when θe = θe,max , the overflow delay achieves its minimum of zero, i.e.g∗ (θ e,max ) = 0. Given that the overflow delay cannot be negative, its associated overflow delay must be either zero or greater, i.e.,g∗(θ e ) ≥ 0 for θ e, min ≤ θ e ≤ θ e, max . On the other hand, the result for Case 1 above indicates that g∗(θ e ) is a piecewise constant function of θ e in each feasibility range. These two results require that at some or all of the breakpoints, there is a sudden drop in the value of g∗(θ e ). The resultant g∗(θ e ) is then a function composed of a piecewise constant segment along each feasibility range, and a sudden drop at each breakpoint. 

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

179

Proposition 4. The phase-2 cost Q(θ e ) is a discontinuous piecewise increasing function with respect to θ e , composing of piecewise increasing segments and breakpoints, segmented in the same way as g∗e (θe ). The minimum Q(θ e ) is achieved at the drop points, which can be identified by the feasibility ranges.

Proof. Case 1: in the feasibility ranges ∗ ∗ Based on the definition of Qe expressed on the RHS of (29), according to Lemmas 2 and 3, i.e. ∂ t∂θ(θe ) = γ and ∂ g∂θ(θe ) = 0, e e ∗ θ )+g∗ (θ )) e e we obtain: ∂ Q∂θ(θe ) = ∂ (θe +t (∂θ = 1 + γ . It is a piecewise linear segment within each feasibility range and an increasing e

e

function if 1 + γ ≥ 0 and decreasing function otherwise. According to (54)–(56), if Zd > 0, we must have γ = − α6 . It d

c

d

indicates that if flexible services are deployed for an OD d, the Lagrange multiplier can be calculated as γ = − α6 . Since the d

c

d

unit flexible services cost cd6 (operating and passenger cost combined) is set to be larger than the passenger travel time cost

α d on RTL, we have cd6 ≥ α d ≥ 0. Hence we have setting in this paper.

∂ Q ( θe ) αd ∂θe = 1 + γ = 1 − c6 ≥ 0. It is an increasing function according to the d

Case 2: at the breakpoints g∗(θ e ) is discontinuous whereas t∗(θ e ) is continuous; hence Q (θe ) = θe + t ∗ (θe ) + g∗ (θe ) is also discontinuous at any breakpoint. Combining the two cases, the phase-2 cost function Q(θ e ) inherits the properties of g∗(θ e ), i.e. piecewise linear and non-convex. The minimum Q(θ e ) must be achieved at the drop points, which can be identified by the feasibility ranges. 

Notation Table.

Parameters N A G(N, A ) (i, j ) ∈ A li j r∈R Rmax fmin , fmax ir ( iu , iv ) ( i u , ju ) S = {Sr } T = {Tr } Sr , Tr (Sr , ir ) (ir , Tr ) N =N∪S∪T A d∈D dmax Bd B¯ d bd

σd d d

n ci1j ci2j c3 ci4j c5 cd6

τ

C M m e∈E pe

node set link set transportation network feasible link from node i to j link distance rapid transit line set the maximum number of lines allowed lower and upper bounds of the line frequency sub nodes of station i: station gate sub node if r = 0, and line sub node if r = 1, . . . Rmax in station link from sub node iu to iv out station link from station i to j on transit line u dummy origin set dummy destination set dummy origin and dummy destination nodes of line r Dummy link connecting dummy origin to line sub node Dummy link connecting line sub node to dummy destination the set of all nodes including the dummy origin and destination nodes the set of all links including the dummy links demand set the number of OD pairs random demand of OD pair d demand covered by RTL in phase-1 for OD pair d demand expectation of OD pair d demand standard deviation of OD pair d cumulative distribution function (CDF) for demand of OD pair d redefined uncertainty set of OD pair d with d = [bd − n ∗ σ d , bd + n ∗ σ d ] breadth of uncertainty set d transit operating cost of link i j transit construction cost of link i j unit station construction cost passenger travel cost on link i j passenger value of time combined flexible services cost of OD pair d consisting of the company operating cost and passenger travel cost transfer penalty transit unit capacity an extremely large positive number number of demand scenarios demand scenario index probability of scenario e happens (continued on next page)

180

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

Decision variables Y = {Yirj } f = { fr } y = {yrij } Y = {Y i j } W = {Wi } Z = {Z d } X = {Xidu jv } d Xboard d Xalight d Xtrans f er d Xin− vehicle ρ = {ρ d }

ρ∗

Evaluation parameters J (ρ ) Qe Q¯ (ρ )

φ (ρ ) βir,e j −βir,e j β = {βirj } d δ ρd fd

fd k kmax

πk λk ,ς θe θe,min ,θe,max [θe , θe ] θe,opt te (θe ) ge (θe ) K d = {κ} hdκ

ηidj,κ αd ω α, β , γ

RTL alignment, Yirj is 1 if link (i, j ) is on line r and 0 otherwise. frequency of line r frequency of link i j on line r RTL alignment Y i j = 1 if link i j is covered by any line r ∈ R station selection, Wi = 1 if any line r ∈ R traverses node i flexible services deployment for OD pair d passenger flow set on RTL for OD pair d through link (iu , jv ) boarding passenger flow subset for i = j, u = 0, v ∈ R in X = {Xidu jv } alighting passenger flow subset for i = j, u ∈ R, v = 0 in X = {Xidu jv } transfer passenger flow subset for i = j, u, v ∈ R, u = v in X = {Xidu jv } in-vehicle passenger flow subset for i = j, u = v ∈ R in X = {Xidu jv } service reliability optimal service reliability RTL construction and operating cost phase-2 cost under scenario e expected phase-2 cost as a function of ρ overall system cost as a function of ρ dual solution corresponding to constraint (31) in the phase-2 problem for link i j on line r of scenario e overflow delay  set of dual solution expectations with βirj = e∈E pe βir,e j a small perturbation in robustness level trial value of service reliability with a small change δ d in ρ d frequency perturbation due to a small change δ d in ρ d the set of frequency difference iteration number the maximum iteration number step size in iteration k parameters in step size calculation in iteration k total flexible services cost for scenario e lower bound and upper bound of flexible service cost feasibility range of θe optimal flexible services cost passenger travel time cost as a function of θe overflow delay cost as a function of θe passenger path set on RTL of OD pair d passenger path flow on path κ of OD pair d Link-path incidence, ηidj,κ = 1 if path κ of OD pair d traverses link i j, 0 otherwise path travel time for OD pair d step size in searching θe,opt dual optimal solution set of P3

References An, K., Lo, H.K., 2014. Service reliability based transit network design with stochastic demand. Transportation Research Record 2467, 101–109. An, K., Lo, H.K., 2014. Ferry service network design with stochastic demand under user equilibrium flows. Transportation Research Part B 66, 70–89. An, K., Lo, H.K., 2015. Robust transit network design with stochastic demand considering development density. Transportation Research Part B 81 (3), 737–754. Bar-Gera, H., 2009. Transportation Network Test Problems. http://www.bgu.ac.il/∼bargera/tntp/ (accessed on Jan 15, 2015) Bell, M.G.H., 1995. Stochastic user equilibrium assignment in networks with queues. Transportation Research Part B 29 (2), 125–137. Ben-Tal, A., Chung, B.D., Mandala, S.R., Yao, T., 2011. Robust optimization for emergency logistics planning: risk mitigation in humanitarian relief supply chains. Transportation Research Part B 45 (8), 1177–1189. Ben-Tal, A., Nemirovski, A., 1999. Robust solutions of uncertain linear programs. Operations Research Letters 25 (1), 1–13. Bertsimas, D., Goyal, V., 2010. On the power of robust solutions in two-stage stochastic and adaptive optimization problems. Mathematics of Operations Research 35 (2), 284–305. Bertsimas, D., Sim, M., 2004. The price of robustness. Operations Research 52 (1), 35–53. Bertsimas, D., Tsitsiklis, J., 1997. Introduction to Linear Optimization (Athena Scientific Series in Optimization and Neural Computation, 6), 1st ed. Athena Scientific ISSN/ISBN: 1886529191. Birge, J.R., Louveaux, F.V., 1988. A multicut algorithm for two-stage stochastic linear programs. European Journal of Operational Research 34 (3), 384–392. Bruno, G., Ghiani, G., Improta, G., 1998. A multi-modal approach to the location of a rapid transit line. European Journal of Operational Research 104 (2), 321–332. Chen, X., Sim, M., Sun, P., 2007. A robust optimization perspective on stochastic programming. Operations Research 55 (6), 1058–1071. Chung, B.D., Yao, T., Friesz, T.L., Liu, H., 2012. Dynamic congestion pricing with demand uncertainty: a robust optimization approach. Transportation Research Part B 46 (10), 1504–1518. De-Los-Santos, A., Laporte, G., Mesa, J.A., Perea, F., 2012. Evaluating passenger robustness in a rail transit network. Transportation Research Part C 20 (1), 34–46. Fan, W., Machemehl, R., 2006. optimal transit route network design problem with variable transit demand: genetic algorithm approach. Journal of Transportation Engineering 132 (1), 40–51. Gabrel, V., Murat, C., Thiele, A., 2014. Recent advances in robust optimization: an overview. European Journal of Operational Research 235 (3), 471–483. Guan, J.F., Yang, H., Wirasinghe, S.C., 2006. Simultaneous optimization of transit line configuration and passenger line assignment. Transportation Research Part B 40 (10), 885–902. Guihaire, V., Hao, J., 2008. Transit network design and scheduling: a global review. Transportation Research Part A 42 (10), 1251–1273.

K. An, H.K. Lo / Transportation Research Part B 84 (2016) 157–181

181

Hensher, D.A., Li, Z., 2012. Ridership drivers of bus rapid transit systems. Transportation 39 (6), 1209–1221. Lam, W.H.K., Gao, Z.Y., Chan, K.S., Yang, H., 1999. A stochastic user equilibrium assignment model for congested transit networks. Transportation Research Part B 33 (5), 351–368. Laporte, G., Marin, A, Mesa, J., Ortega, F., 2007. An integrated methodology for the rapid transit network design problem. Algorithmic Methods for Railway Optimization. Springer, Berlin Heidelberg, pp. 187–199. Laporte, G., Mesa, J.A., Perea, F., 2010. A game theoretic framework for the robust railway transit network design problem. Transportation Research Part B 44 (4), 447–459. Laporte, G., Mesa, J., Ortega, F., Sevillano, I., 2005. Maximizing trip coverage in the location of a single rapid transit alignment. Annals of Operations Research 136 (1), 49–63. List, G.F., Wood, B., Nozick, L.K., Turnquist, M.A., Jones, D.A., Kjeldgaard, E.A., Lawton, C.R., 2003. Robust optimization for fleet planning under uncertainty. Transportation Research Part E 39 (3), 209–227. Lium, A., Crainic, T.G., 2009. A study of demand stochasticity in service network design. Transportation Science 43 (2), 144–157. Lo, H.K., An, K., Lin, W., 2013. Ferry service network design under demand uncertainty. Transportation Research Part E 59, 48–70. Lou, Y., Yin, Y., Lawphongpanich, S., 2009. Robust approach to discrete network designs with demand uncertainty. Transportation Research Record: Journal of the Transportation Research Board 2090, 86–94. Marin, A., Jaramillo, P., 2008. Urban rapid transit network capacity expansion. European Journal of Operational Research 191 (1), 45–60. Meng, Q., Lam, W.H.K., Yang, L., 2008. General stochastic user equilibrium traffic assignment problem with link capacity constraints. Journal of Advanced Transportation 42 (4), 429–465. Patriksson, M., 2008. On the applicability and solution of bilevel optimization models in transportation science: a study on the existence, stability and computation of optimal solutions to stochastic mathematical programs with equilibrium constraints. Transportation Research Part B 42 (10), 843–860. Perea, F., Mesa, J.A., Laporte, G., 2014. Adding a new station and a road link to a road–rail network in the presence of modal competition. Transportation Research Part B 68, 1–16. Samanta, S., Jha, M.K., 2011. Modeling a rail transit alignment considering different objectives. Transportation Research Part A 45 (1), 31–45. Slyke, R.M.V., Wets, R., 1969. L-shaped linear programs with applications to optimal control and stochastic programming. SIAM Journal on Applied Mathematics 17 (4), 638–663. Soyster, A.L., 1973. Technical note—convex programming with set-inclusive constraints and applications to inexact linear programming. Operations Research 21 (5), 1154–1157. Walteros, J.L., Medaglia, A.L., Riaño, G., 2015. Hybrid algorithm for route design on bus rapid transit systems. Transportation Science 49 (1), 66–84. Wan, Q., Lo, H., 2009. Congested multimodal transit network design. Public Transport 1 (3), 233–251. Wan, Q., Lo, H., 2003. A mixed integer formulation for multiple-route transit network design. Journal of Mathematical Modelling and Algorithms 2 (4), 299–308. Wang, S., Meng, Q., 2012. Robust schedule design for liner shipping services. Transportation Research Part E 48 (6), 1093–1106. Yan, Y., Meng, Q., Wang, S., Guo, X., 2012. Robust optimization model of schedule design for a fixed bus route. Transportation Research Part C 25, 113–121. Yu, B., Yang, Z., Jin, P., Wu, S., Yao, B., 2012. Transit route network design-maximizing direct and transfer demand density. Transportation Research Part C 22, 58–75.