An integrated model of energy-efficient timetabling of the urban rail transit system with multiple interconnected lines

Urban rail transit (URT) has been considered an effective means of addressing urban congestion problems in metropolises. The operations of a URT system involve high energy consumption and its trade-off with passenger travel times. Existing energy-efficient timetabling studies have predominately focused on single URT lines and thus are incapable of accurately modeling the energy consumption in a URT network with transfer opportunities and synchronization between the URT lines. To extend the energy-efficient timetabling from one single line to multiple interlinked lines, we propose a bi-level model incorporating the operator’s decision on a regular timetable and passengers’ path choice in a URT network. The objective of energy consumption and timetable constraints of the upper level are linearized and formulated as mixed-integer linear programming. The lower level captures the user equilibrium based path choice behavior responding to the timetable. We develop a heuristic algorithm for the bi-level model that produces near-optimal timetable solutions in a relaxation process. The proposed model and solution algorithm are validated in the URT network of Xi’an (China). It is found that the energy consumption is considerably reduced, compared with using the current timetable, at the expense of an acceptable increase in the average travel time.


Introduction
Urban rail transit (URT), due to the low fare, high capacity, punctuality, and environmental friendliness, has attracted much attention in urban mobility in large cities . However, even though the URT is much more energy-efficient than other main transport modes for serving the same population (González-Gil et al., 2013, the URT faces a challenge in energy-saving. For instance, it is estimated that the Beijing URT will approximately use two billion kWh in 2020 (the same amount may supply a city with a population of 500 thousand for one year) and be the largest consumer of all single organizations in Beijing, which cause concerns to the local government and environment agencies (Lv et al., 2019). In the existing research, there are two research lines dedicated to energy efficiency in a URT system , namely, optimal train control and energy-efficient timetabling. Focusing on a single train on a track between two stations, the optimal train control aims at finding the optimal driving strategies in terms of speed profiles to minimize energy consumption (Howlett, 1996(Howlett, , 2000Howlett and Pudney, 1998;Albrecht et al., 2013;Scheepmaker et al., 2017). On the other hand, energy-efficient timetabling (Li and Lo, 2014a, b;Yang et al., 2015Yang et al., , 2016Yang et al., , 2017Yin et al., 2016;Goverde, 2016, 2017) typically concerns the running time allocations of a fleet of trains on a URT line. As energy-efficient timetabling is the focus of attention, we review the most relevant studies in two research lines, i.e., (1) optimal train control for reducing energy consumption and (2) energy-efficient timetabling for single metro lines.
Optimal train control is a traditional problem based on the optimal control theory, in particular on the Pontryagin's Maximum Principle (PMP) (Albrecht et al., 2016a,b). The mathematical problem was conducted and applied under the control theory in the early years (Howlett and Pudney, 1995;Khmelnitsky, 2000;Liu and Golovitcher, 2003;Lai et al., 2020). The train control and speed profile include three sub-processes, namely, maximum acceleration, coasting, and maximum braking (Howlett and Pudney, 1998). The switching points of these sub-processes are at the core of train control. Howlett et al. (2009) showed that the optimal switching points for a steep section can be found by minimizing an intrinsic local energy function. Yang et al. (2018) analyzed the switching points and converted the relationship between running time and energy consumption into a strictly convex quadratic programming problem. To assist train drivers with optimal controls in real-time, some train on-board systems incorporated specific speed advice and train delays in operations (Panou et al., 2013). With the mass data generated in real operations and development of artificial intelligence, Huang et al. (2019) and Yin et al. (2020) applied machine learning methods, for example, random forest regression (RFR), support vector machine regression (SVR), and deep neural network (DNN), to optimize speed profiles.
Energy-efficient timetabling has recently been a hotspot topic (Li and Lo, 2014a, b;Zhao et al., 2015;Gupta et al., 2016;Ye and Liu, 2016;Yin et al., 2016Yin et al., , 2017Canca and Zarzo, 2017;Yang et al., 2019;Mo et al., 2019a,b;Liu et al., 2020;Qu et al., 2020) in the study of energy efficiency of URT systems. For example, Li and Lo (2014a) considered the speed control and headway of a timetable by synchronizing train acceleration and braking to maximize the utilization of regenerative energy. It was later found by Li and Lo (2014b) that adjusting the cycle time of the URT line could further reduce energy consumption. These two studies formulated the problems in a non-linear system and applied a genetic algorithm (GA) to find the solutions. Considering uncertain dwell times, Yang et al. (2017) proposed an ε-constraint method in the GA framework to find Pareto optimal solutions to minimize the total energy consumption and passenger travel time. Utilizing passenger smart-card data as input for passenger demand, Yang et al. (2019) formulated a convex quadratic programming problem and developed an optimization-based approach to generate a timetable with stop-skipping patterns. By transferring arrival and departure times to time window constraints and relaxing the given timetable, Wang and Goverde (2019) proposed a multi-train trajectory optimization method with a base objective of minimizing multi-train energy consumption and an additional objective of eliminating conflicts between trains. Mo et al. (2019a) presented a bi-objective model to minimize energy consumption and passenger waiting time simultaneously, considering different train capacity constraints. A modified tabu search algorithm with a prior enumeration process was used to find near-optimal solutions. A few factors in a URT network affect the total energy consumption, e.g., rolling stocks and line lengths. Specifically, different lines in a URT network usually operate with different rolling stocks (Zhong et al., 2019), which affects train deadhead mileage, passenger flow, and energy consumption on the sections (Zhong et al., 2020). Mo et al. (2019b) integrated the optimal train timetable and rolling stock plan that optimized the braketraction overlapping time at stations. Based on energy-regenerative technologies, Yang et al. (2020) formulated an optimization model considering energy allocation and passenger assignment to maximize regenerative energy utilization and minimize passenger travel time. A solution algorithm based on the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) was developed to find the Pareto frontier. All of the above studies reported significant energy reductions of the newly generated timetables over the used timetables.
It is commonly known that the timetabling for a single URT line is an NP-hard problem (Cai and Goh, 1994). The timetabling problem becomes harder when energy-efficiency is incorporated. As summarized in Table 1, with the assumption of fixed passenger arrival rates and the absence of path choices, a few studies (e.g., Yin et al., 2017) suggested mixed-integer linear programming (MILP) formulations for single URT lines of less than 20 stations. With a similar assumption, some studies (e.g., Yang et al., 2020) considered non-linear formations and applied metaheuristic solution algorithms (e.g., GA and PSO). The energy-efficiency timetabling problem is even harder when considering passenger path choices in a URT network of multiple interconnected lines, which is probably the reason why no research of this kind has been done.  Li and Lo, (2014a, b) S No NP GA Yang et al. (2015Yang et al. ( , 2016Yang et al. ( , 2017 S No NP GA Gupta et al. (2016) S No NP N/A Wang and Goverde, (2016 S Despite increasing research entries on energy-efficient timetabling, to the best of our knowledge, existing studies have been unexceptionally limited to single URT lines. This treatment is valid in a URT system with only one URT line or independent URT lines, which however does not hold in most modern URT systems. Although the power supply for each URT line is separated, the URT lines are associated with each other by passenger transfers and timetable synchronization among the URT lines. The modeling limitations cannot capture the accurate passenger allocation in the URT network and tend to produce unwanted chain effects on the calculation of energy consumption. Many studies have addressed the traditional timetabling problems of a URT network with multiple lines (e.g., Wong et al., 2008;Guo et al., 2018;Shang et al., 2018;Sun et al., 2018;Chen et al., 2019;Yao et al., 2019;Guo et al., 2020;Boroun et al., 2020;Huang et al., 2021), but none of them takes energy-efficiency into consideration.
Therefore, this study aims to extend the energy-efficient timetabling from a single line to multiple lines in a URT network. We propose a bi-level model framework considering timetabling and passenger path choices (including transfer choices) as an interactive process. The upper level model determines the optimal timetable and speed profiles under the given passenger path choices. The lower level model concerns passenger path choices by the receiving timetable from the upper level. Technically, we formulate the upper level optimization problem as mixed-integer linear programming (MILP). The lower level adopts a passenger assignment in a space-time network and is assumed that passengers adapt path choices until a user equilibrium (UE) state is reached.
To highlight the contributions of this study, we summarize the differences between this study and a few relevant studies in terms of modeling aspects in Table 1. First of all, our study is the first one extending the energy-efficient timetabling from a single line to multiple lines in a URT network. Second, although passenger assignment has recently been considered in energy-efficient timetabling (e.g., Su et al., 2019;Yang et al., 2020), these studies applied a one-off allocation based on simple heuristics of proportional assignment. This study considers more realistic UE-based feedback of passengers to the timetable. Third, most studies formulated one level of (non-)linear programming and applied (meta-)heuristic methods (e.g., GA and PSO), which might result in non-justifiable and non-reproducible solutions. In this study, we propose a bi-level model framework and develop a domain-knowledge based heuristic algorithm to obtain near-optimal solutions. The remainder of this paper is organized as follows. Section 2 provides the problem description and modeling assumptions. Section 3 entails the modeling of synchronous URT lines and their interactions with passenger path choices in a bi-level model framework. Section 4 discusses the solution algorithm to the bi-level model. Section 5 presents a comprehensive case study in the URT system in the city of Xi'an (China) to verify the model and algorithm. Finally, Section 6 concludes the main contributions and provides suggestions for future research.

Problem description
The energy-efficient timetabling determines the elements of a timetable (e.g., arrival time, running time, departure time for each platform) with the minimum energy consumption while completing the tasks of transporting passengers. The passengers in the URT network may have multiple alternative paths. Therefore, passenger path choices should be taken into consideration. In Fig. 1, we use a typical and simplest URT network with two crossed lines to illustrate the timetabling problem. Suppose there are N l and N l ' stations, 2N l and 2N l ' platforms, and 2N l − 2 and 2N l ' − 2 tracks (one track between two neighboring platforms) on line l and l ' , respectively. The platforms and tracks are arranged and numbered in order. A train departing from platform 1 to platform N l is in the up-direction, while departing from the platform N l +1 to platform 2N l is in the down-direction. There are four platforms on the transfer station, i.e., platform 4 and 2N l ' − 3 of line l ' , platform 3 and 2N l ' − 2 of line l. Based on this layout, key components for energy-efficient timetabling are represented in Fig. 2. A passenger starts from its origin platform and goes through a series of activities (waiting, boarding, and alighting) to reach the destination platform.
For illustration, as shown in Fig. 3(a), the running process on each track has a set of alternative discrete speed profiles and each profile has three main stages (accelerating, coasting, and decelerating) with speed limitations for different geographical conditions. Commonly, the speeds and running times are dynamic according to the operation level of the on-board control systems called ATO (Automatic Train Operation). The drivers intervene and control the trains only in case of emergencies (Yin et al., 2017;Mo et al., 2019a). Based on the train control theory, the energy consumption on track t is a non-linear function of the loading weight (m l,t ) and running time n l,t , denoted by e ( m l,t , n l,t ) . For each track, one level of on-board operation corresponds to one speed profile and thus one specific running time. A higher operation level usually corresponds to a lower speed. Denote the running time at operation level g on track t of line l by n l,t,g . For simplicity, n l,t,g is simplified as n g given a track. To illustrate the relationships, suppose that operation level g is associated with running times n g . Given train loading weight level m, the energy consumption e gm (for the sake of notational   simplicity) has the following features. As shown in Fig. 3(b), for the same loading, a higher speed profile needs more energy consumption for acceleration and coasting (the main stages of energy consumption). Therefore, running times and energy consumption satisfy n 1 ≤ n 2 ≤ n 3 and e 11 ≥ e 21 ≥ e 31 respectively. For the same operation level, the heavier the train is, the more energy consumption is involved (i.e., e 12 ≥ e 11 ; e 22 ≥ e 21 ; e 32 ≥ e 31 ).
Overall, operation level g and train loading m l,t determine the running time and energy consumption given operation conditions (train type and ATO system). The train loading, m l,t , depends on the passenger volume on the train. Therefore, there are two types of decision variables, i.e., timetable and passenger distribution on the network, in the energy-efficient timetabling problem for multiple lines. The notations are defined below.

Assumptions
The following assumptions are made to facilitate the development of the model framework.
Assumption 1. The trains are scheduled periodically (e.g., during peak hours) and all passengers can board the train in a planning period. Periodic timetables are widely adopted by URT operators using the space-time network modeling method (Shang et al., 2018.

Assumption 2.
In the URT network, the transfer time consists of walking time and waiting time, where the walking time is fixed given the transfer station and the waiting time is related to the headway of the URT line.

Assumption 3.
To accommodate all passengers, the passenger arrival rate determines the headway for a periodic timetable. The time-dependent passenger arrival rate at a platform is fixed as the average during the planning period.
Assumption 4. Given a timetable, the passengers in the URT network adapt path choice to reach a UE state by long-term day-to-day adjustments (Bie and Lo, 2010;Djavadian and Chow, 2017).
Assumption 5. The URT operator pursues the minimum energy consumption, while the passengers seek the least generalized travel cost.

Notations
The following parameters, sets, intermediate variables, and decision variables are used in the model. Parameters: Elements and sets:

Model
The primary difference between single-line and multiple-lines is whether the transfer is allowed. Passenger path choices are the most salient differences between a single line and multiple lines in a URT network. Although the power supply for each URT line is separated, the URT lines are associated with each other by passenger transfers and synchronization between the URT lines. For a URT network with transfer opportunities, passenger transfer and path choice should be incorporated in the energy-efficient timetabling. Passengers' path choice and the operator's timetable have a typical game relationship (Assumption 5). Therefore, we propose a bi-level model for the energy-efficient timetabling in a URT network to model the interactions. The energy-efficient timetabling model at the upper level is formulated as mixed-integer linear programming. The lower level is specified as a UE-based passenger assignment in the URT network. Fig. 4 shows the bi-level model framework, in which the upper level is an energy-efficient timetabling model that needs the input of passenger path choices, and the lower level is the passenger assignment responding to a timetable. The bi-level optimization framework divides the decision variables into two parts, i.e., timetable and passenger flow patterns in the URT network. Since it is not straightforward to formulate the lower-level user equilibrium conditions as constraints in the upper-level model, the timetabling (up-level) and passenger flow assignment (lower-level) are carried out in a heuristic relaxation process to seek quality timetable solutions.

Energy-efficient timetabling model (Upper level)
Focusing on transfer between multiple lines, we optimize a periodic timetable with a fixed headway of each line at the upper level. In the following two subsections, we first propose the energy-efficient timetabling in non-linear formulations to minimize the total energy consumption. Then, the model is transformed into MILP by linearizing the non-linear objective and constraints.

Energy-efficient timetabling model in non-linear formulation
The timetable-related constraints and objective of energy consumption are presented as follows. As the designed timetable is periodic and parallel, the time schedules of the first train are repeated cyclically by the ensuing trains on the same URT line. For track t on line l, ∀t ∈ T l , ∀l ∈ L, the definitional relationship among departure time (d l,t ), running time (n l,t ) and arrival time (a l,t+1 ) is formulated as Eq. (1). For the dwell time on platform p (∀p ∈ P l ) on line l, the relationship between departure time (d l,p ) and arrival time (a l,p ) is formulated as Eq. (2). To allow enough time for alighting and boarding and avoid waiting long at the platform, there is a time range, The minimum headway is determined by the signaling system for safety considerations and the maximum value is usually related to the service level. The departure frequency f l determines the train fleet size during the planning period and is associated with the headway by Eq. (5). To ensure frequency f l is an integer, h l is specified as a divisor of 3600 (considering one hour as the time frame). As Canca and Zarzo, (2017) suggested, there are alternative headway values, for instance, 120, 180, 240, 300, 360, 600, 720, 900, 1200, and 1800 in seconds. The cycle time Φ l is the time expense by a train completing a loop in the bi-directional line, given as Eq. (6), where η l denotes the turning-around time of line l.
The passenger volume on track t can be computed by the accumulations on the track, which is illustrated in Fig. 5. Starting from the first platform in one direction, v + l,p and v − l,p denote the numbers of boarding and alighting passengers at platform p of line l respectively. As shown in Eq. (7), the passenger volume on track t is The dwell times need to be long enough to complete the alighting and boarding processes. As shown in Fig. 6, the alighting and boarding processes take time . The dwell time should be no less than the alighting and boarding time, i.e., b l,p ≥ ( y − l,p +y + l,p ) , rearranged as (12). Since the opening and closing times of the train doors are short, they are not considered in the dwelling process.
The headway is associated with the fleet size φ l , expressed as Eq. (13) as Eq. (14). We can see that if cycle time Φ l ≥ 3600s, then φ l ≥ f l . If Φ l < 3600s (there may be URT lines with short cycles), then φ l < f l . Economically, the fleet size should not be greater than the maximum number of trains in Eq. (15). There is a corresponding capacity of each line, γ l , for the consideration of limited room and passenger comfort. To accommodate all passenger demands given the capacity, frequency f l should be large enough to satisfy the constraint as Eq. (16). Constraint (16) can be in another form due to f l •h l = 3600 as Eq. (17). Therefore, if the train capacity is fixed after determining the train units, the maximum passenger volume in one hour is 3600•γ l h E . Additionally, with an increase in q max l , the headway should be shorter. Given average passenger weight τ, the passenger load m l,t on track t is formulated as Eq. (18).
The objective of energy-efficient timetabling is the total energy consumption of the URT network. The energy consumption for line l is given as Eq. (19) and the total energy consumption of the URT network is given as Eq. (20).
where e l is the energy consumption of linel and E is the total energy consumption of the URT network. As formulated in Eq. (20), the energy consumptions of the URT lines are independent of each other. However, the linkage between the URT lines are reflected in q l,t of Eq. (10) and m l,t of Eq. (19) during the course of passenger path choice. One derived feature is that Eq. (20) can be solved in parallel given the passenger path choices. This implicit decomposition of energy consumption reduces the complexity of energy-efficient timetabling.

Energy-efficient timetabling model as MILP
To make the energy-efficient timetabling problem tractable, we linearize the non-linear objective function (Eq. (20)) and two constraints (Eqs. (5) and (13)). First, the non-linear objective function is linearized by utilizing discrete speed levels. Then, based on alternative headways, the constraints are linearized by the Big M method and variable transformations.
Denote the basic energy consumption on track t at level g by e 0 l,t,g when the train is empty with weight m 0 l . When the passengers are loaded on the train, the total weight is m 0 l + m l,t . Thus, the energy consumption in a neat form is equal to e 0 l,t,g (the multiplication sign • is omitted for simplicity). The side-effect of this formulation on the energy consumption is negligible because the energy consumption is mainly determined by the mechanical running resistance and the kinetic energy (other parts, such as aerodynamical resistance and the ingoing air volume, are only in very small proportions). Given a set of operation levels on track t of line l, G lt , a set of binary variables θ l,t,g is introduced: θ l,t,g = 1 if level g is selected; otherwise, θ l,t,g = 0. For track t, only one level can be selected, expressed as Eq. (21). Given the selected level, the running time on a track is fixed. Therefore, the running time for track t of line l can be formulated in Eq. (22). ∑ g∈Glt θ l,t,g = 1, ∀l ∈ L, ∀t ∈ T l (21) Therefore, energy consumption e ( m l,t , n l,t ) in Eq. (19) can be replaced by Eq. (23) consists of the product of binary variable θ l,t,g and variable . To linearize Eq. (23), we introduce a new set of variables ρ l,t,g : ρ l,t,g = if θ l,t,g = 1; otherwise, ρ l,t,g = 0. Then, a set of constraints are introduced as where Λ 1 is an upper bound of ρ l,t,g , for instance, max Denote an element of H l by h l,k ( k ≤ |Ω l | ) and an element of Ω l by f l,k . To keep h l,k f l,k = 3600, h l,k is in ascending order if f l,k is given in descending order.
To linearize the constraints, we introduce a series of intermediate binary variables δ l,k . The sum of the binary variables δ l,k is equal to 1 in Eq. (25). The binary variable δ l,k and the element h l,k can be utilized to represent the headway as the sum of their product in Eq. (26). Similarly, the frequency can be expressed in Eq. (27). Recall that constraints (5) and (13) are non-linear. After introducing δ l,k , H l , and Ω l , constraint (5) is linearized. Constraint (13) can be rewritten as Eq. (28).
To linearize Eq. (28), a set of intermediate variables α l,k is further introduced: α l,k = φ l if δ l,k = 1; otherwise, α l,k = 0. The linearization is achieved by a set of constraints in Eq. (29), where Λ 2 is an upper bound of α l,k , for example, Λ 2 = maxφ l h l ,∀l ∈ L. The cycle time of line l is given in Eq. (30).
With the new variables, the objective is expressed as Furtherly, we replace (ρ l,t,g δ l,k ) by a new variable λ l,t,g,k , where λ l,t,g,k = ρ l,t,g when δ l,k = 1 and λ l,t,g,k = 0 when δ l,k = 0. The relationship between λ l,t,g,k , ρ l,t,g , and δ l,k are described by constraints (32), where Λ 3 is the upper bound of λ l,t,g,k , for instance, Taken together, the linear formulation of the objective is shown as Eq. (33).
≤ λ l,t,g,k ≤ ρ l,t,g + Λ 3 ( 1 − δ l,k ) λ l,t,g,k ≤ δ l,k Λ 3 λ l,t,g,k ≥ 0 ∀l ∈ L, t ∈ T l , g ∈ G lt , k ≤ |Ω l | With the objective and constraints linearized, the upper level model (UL for short) is transformed into a MILP formulation as As the UL is in mixed-integer linear formulations, we can obtain the exact energy-efficient timetable solution by optimization solvers (e.g., Cplex and Gurobi) given the passenger path choices. It is noteworthy that the linearization of the constraints by the Big M method and variable transformation does not modify the space of feasible solutions. In the UL, the total number of variables is the sum of those in each line. The detailed numbers of variables and constraints are listed in Table 2 to reflect the true complexity of the model. The total number of variables is 4|L| + 3 where operator |•| gives the cardinality of a set. The total number of constraints is 7|L| + 2 ∑ l |P l | + 8 The numbers of variables and constraints for line l are 3|P l | +2|T l | +|Ω l | + ∑ t |G lt ||Ω l |, ∀t ∈ T l , respectively. Due to the implicit decomposition, the URT line with the most variables and constraints dictates the complexity of the UL, which is under the capacity of the optimization solvers.

Passenger assignment (Lower level)
Passenger path choice is the key to URT operations. Much research has been done on passenger path choice in traffic networks. The majority of studies revolve around the UE-based traffic assignment (Wardrop, 1952;Nguyen and Pallottino, 1988;Spiess and Florian, 1989;Wu et al., 1994;Nuzzolo et al., 2001;Gentile et al., 2005;Jiang and Szeto, 2016;Shang et al., 2019,). However, little research has been thus far dedicated to passenger assignment for energy-efficient timetabling in a URT network. We formulate a passenger assignment model in a one-time instance of the space-time URT network as the lower level (LL for short) of the bi-level model.
We first construct the space-time URT network responding to a timetable generated from the UL model. In the space-time URT network, crowding, waiting, and transfers should be considered for defining passenger travel cost (disutility) on the paths. Given a timetable, the passenger flow patterns influence crowding on different running arcs. In line with the literature of path choice in traffic networks, we assume that the passengers adjust path choices to reach a user equilibrium (UE) state, at which all selected paths have an equal and minimum cost (Assumption 4). The space-time network consists of five types of arcs (virtual arc, waiting arc, running arc, transfer arc, and arrival arc) and two types of nodes (physical node for the platform and virtual node for the station). One station is extended into one virtual node and two real nodes corresponding to the two platforms. The virtual arc connects the platforms of the same station. The waiting arc represents the waiting stage before boarding at the origin station of a trip. Running arcs represent physical movements over the tracks. A transfer arc bridges the platforms of the feeder line and the platforms of the connecting line at the transfer station. The arrival arc refers to the final stage of a trip leaving the destination platform. One temporal instance of the space-time network is illustrated in Fig. 7. For passengers departing from origin node o to destination node d, they first traverse the waiting arc to take the train on the running arc; after that, they arrive at the transfer station and go through the transfer arc to the connecting line. Next, they take running arcs again and finally reach the destination. Since we focus on a periodic timetable (Assumption 1), the space-time network has a common backbone structure. The generalized costs of the different arcs are defined below.
The cost of the virtual arc i, c i (i ∈ I v ), is considered zero in Eq. (34). Similarly, the cost of the arrival arc c i (i ∈ I a ) is also considered zero.  (1) where I v and I a are the sets of virtual arcs and arrival arcs respectively. According to Assumption 4, the time-dependent passenger arrival rate at a platform is fixed as the average during the planning period. The average waiting time at a platform is considered half of the headway (Ceder and Tal, 2001). Thus, the cost of a waiting arc is simply set as where κ w is the cost coefficient of waiting time and I w is the set of waiting arcs. In the literature of route choice in a public transport system, crowding is often considered an important factor affecting passenger route choice Fu and Lam, 2018). In this study, the passenger crowding effect is considered on the running arcs of the URT network. The crowding is directly determined by the passenger volume. Note that the crowding in URT is different from the crowding on the road surfaces. Crowding on the roads is out-the-vehicle which causes losses of time, but crowding in URT is in-vehicle which causes discomfort other than losses of time. Therefore, the cost of a running arc is determined by the running time of each track and the crowding effect , formulated in a BPR-like (Bureau of Public Road) form as where κ r is the crowding coefficient of in-vehicle time and I r is the set of running arcs. According to Assumption 2 and Eq. (36), the cost of a transfer arc is given as  where κ t is cost coefficient of transfer time, w i is the fixed walking time for transfer arc i at a specific station, h i is the headway of the connecting line on transfer arc i, and I t is the set of transfer arcs. Based on the generalized arc costs in the space-time URT network, the generalized path cost for an OD pair is formulated in an additive form as Eq. (38). The UE condition is given in Eq. (39).
Eq. (38) is the additive costs of a path, where u ijr is a 0-1 variable, u ijr = 1 if the arc i belongs to path r of OD pair j, otherwise, u ijr = 0. Eq. (39) describes the equilibrium condition, i.e., the used paths for an OD have equal and minimum cost and unused paths have costs higher than or equal to the minimum cost. μ j is the minimum cost of OD j at the UE state, z r j is the cost of path r that belongs to OD j, and q r j is the passenger flow assigned to path r of OD j. Eq. (40) is the flow conservation on paths, where Q j is the travel demand of OD j. Eq. (41) gives the passenger volume of arc i. To get the passenger volume q l,t on the track, β l,t i (a 0-1 variable) is introduced in Eq. (42) to represent the passenger volume of arc i, where β l,t i = 1 if arc i corresponds to track t on line l, otherwise, β l,t i = 0. In any time instance of the space-time URT network, all expanded arc costs are static and thus the path costs are continuous and non-decreasing with path inflows according to Eqs. (34)-(38). Therefore, given a periodic timetable, a UE state exists after long-term adaptions of path choice (including transfer choices). Existing algorithms such as the Method of Successive Average (MSA) and Routing-Swapping Method are sufficient to find a UE state. The results of the passenger choices at the UE state will be fed to the UL model through Eq. (42) in the bi-level model.
Overall, the proposed optimization framework for generating a periodic energy-efficient timetable relies on assumptions that all passengers can board the train (Assumption 1) and the time-dependent passenger arrival rates at a platform are fixed (Assumption 3) during this period. These assumptions hold in most cases when the passenger demands are evenly distributed in the temporal dimension and not extremely high in the planning period. The model has limitations with uneven and large passenger demands. In that sense, the model framework can be extended by incorporating aperiodic timetables and passenger delays, which however substantially increases the solution space and the difficulty to find a quality solution. As the first endeavor to extend the energy-efficient timetabling from a single line to multiple lines, the bi-level optimization framework divides the decision variables into two parts (i.e., timetable and passenger flow pattens) to cut down the complexity. As it is not straightforward to incorporate the lower-level UE condition as a constraint in the upper-level, we have to resort to a heuristic relaxation process, which cannot guarantee obtaining an optimal solution.

Solution algorithm
Bi-level programming is widely applied in transportation research considering two-player interactions (Bracken and McGill, 1973;Yang and Huang, 2004;Gao et al., 2005;Han et al., 2015;Yu et al., 2015;Rashidi et al., 2016;Li and Wan, 2019;Li and Liao, 2020). The bi-level programming for network design problems is proved as an NP-hard (NP: non-deterministic polynomial-time) problem. Several algorithms have been proposed to solve the bi-level model, for example, meta-heuristic algorithms and Lagrangian relaxationbased algorithms (Cantarella et al., 2006;Szeto and Jiang, 2014;Liu and Zhou, 2016). As remarked in Section 3.1, given the passenger path choices, we can directly apply an existing solver to obtain the exact solutions. Since the lower level is a static user equilibrium in the space-time network, the method of successive average (MSA) as a proven method is applied. There are usually four steps of the MSA (see the pseudo-code below), including initialization, path generation, traffic impedance updating, and passenger flow assignment (Mounce and Carey, 2015). Although the MSA may get a local optimum, the outer iterations of the upper and lower level models can be utilized to evaluate the quality of a feasible solution.
The pseudo-code of the MSA: Step 1: Initialization Set passenger flow vector x as 0; calculate traffic impedances for all arcs by Eqs. (34)-(38); set iteration counter s = 0.

Step 2: Path generation
Assign passenger demands to the feasible paths generated by the shortest path algorithm with the minimum path impedances; obtain path and arc flows.
Step 3: Update arc impedances s = s + 1; calculate the traffic impedances of all arcs with passenger flow x s .

Step 4: Update the passenger flow assignment
Repeat Step 2; obtain new passenger flow y s and update x s+1 = x s + (y s − x s )/s.

Step 5: Termination condition
If s ≥ A or ‖x s+1 − x s ‖/x s ≤ ω, then stop, where A and ω are the maximum number of iterations and threshold for convergence respectively; otherwise, s = s +1 and return to Step 3.
For the bi-level model, any feasible timetable and passenger path flow together generate an upper bound. To examine the quality of the upper bound, we propose a domain-knowledge based heuristic method to find a lower bound of the optimal energy consumption. Similar to the MSA algorithm performed in the space-time URT network, the MSA algorithm is carried out in a space-energy network to obtain the lower-bound energy consumption. The space-energy network has the same topology as a one-time instance of the space--time network. Link impedances in the space-energy network refer to energy consumption and the energy-related MSA algorithm is called EMSA for short. One major difference is that only running arcs consume energy and the other four types of arcs do not. To obtain a valid lower bound, we use the profile with the minimum energy consumption of an empty train, e 0,min l,t , on a running arc. For a further relaxation of the lower bound, we only consider a proportion of the train passengers obtained from the EMSA. It is supposed that the passengers on the space-energy network also reach equilibrium with the lowest energy consumption. Likewise, the impedance of running arcs in the EMSA is updated as where q i is the passenger volume on running arc i.
Using the minimum energy consumption on each track, a relaxed constraint for the fleet size, and a similar MSA convergence, the EMSA can find a quality approximate lower bound. Integrating the above modeling components, a solution algorithm to the bi-level model is shown in Fig. 8. Note that, starting from a timetable input, the passenger assignment should be conducted before energy-efficient timetabling. Therefore, the lower level module is swapped to the upside with reference to Fig. 8. The solution algorithm has two stages, in which the "Lower bound" module on the left-hand side aims to find a quality lower bound in Fig. 8(a) and the counterpart on the right-hand side finds a feasible solution for the upper bound in Fig. 8(b). Specifically, at the first stage, we construct the space-energy network S L and obtain the passenger flow X L at an imaginative UE state by EMSA; then, we solve the energy-efficient model Eq. (33) to find the lower bound energy consumption E L and timetable B L . This stage is completed when a converged E L is found. Likewise, at the second stage, we construct the space-time network S U , obtain passenger flow X U by MSA, and find the minimum energy E U with feasible timetable B U . If the relative gap, EU− EL EU , between the upper and lower bounds is not satisfactory, B U is considered the input to the "Upper bound" module to start the next iteration. The solution algorithm is terminated when a satisfactory relative gap is obtained or the maximum number of iterations is reached.

Case study
In this section, we demonstrate the effectiveness of the proposed bi-level model with the URT network in the City of Xi'an (China), which includes 4 lines and 94 stations in service (Fig. 9). Lines 1-4 have 38, 42, 52, and 56 bidirectional platforms respectively. The alternative headways for each line are listed in Table 3. Since a station is extended to three nodes (i.e., two platforms and one virtual node), there are 282 (or 94 × 3) nodes in a one-time instance of the space-time network. We select a URT operating period between 8:00 am and 9:00 am (Γ = 3600s) during the morning peak time. The running time and speed profiles for each section are listed in Table A1 (in the Appendix). The model is solved with a personal computer (8G RAM and Intel Core i7-6700U CPU) and Gurobi + python. The parameters in the lower level UE model are set in the unit of disutility/min as, the cost coefficient of waiting time κ w :1.0, the cost coefficient for crowding of in-vehicle time κ r : 0.1, cost coefficient range of transfer time κ t : [1.3, 2.5] (a feasible range suggested by Chen et al. (2015)). The maximum iterations (A) for MSA and EMSA are set as 500, and ω is set as 10 − 3 . The relative gap threshold, gap 0 , is set as 5% for the bi-level model. To verify the model and algorithm, two cases are presented in an accumulative way. Case 1 shows the effectiveness of the model in the current URT network with a sensitivity test on the key parameter κ t . To highlight the effects of transfer opportunities in the URT network, Case 2.1 supposes disruptions at the transfer stations; in Case 2.2, we add transfer opportunities from and to new URT lines to be opened to study the influence on energy consumption.

Case 1: URT in the current form
Based on the passenger demand on a weekday, the original timetable consumes energy of 94453.98 kWh, calculated by our model. Considering different values of κ t and 90% of the train passengers in the lower bound (Section 4), stable timetable solutions are obtained as shown in Fig. 10. Since all upper bounds are stable after ten iterations, the convergent solutions are shown for the first ten outer iterations, which take around 650 s as the average computation time. Taking κ t = 1.3 for example, a common setup in MSA , the upper bound timetable solution consumes energy of 69537.97 kWh as opposed to the lower bound 66807.877 kWh, resulting in a relative gap of 11.97% (or (69537.97-66807.877)/69537.97). Compared with using the current timetable, the solution Table 3 The alternative headways for each line.      K. Huang et al. timetable reaches a relative energy reduction by 26. 4% (or (94453.98-69537.97)/94453.98). When the κ t is bigger, it may obtain more energy reduction. For example, κ t = 1.5, the relative energy reduction reaches 29.3%. Due to the limited alternatives of train headways (Table 3) and the use of a conservative low bound, the solution space is not large in terms of discrete combinations. We find that the upper bounds are convergent after two iterations and the relative gaps do not reach the pre-specified level with certain specifications of κ t . The reason is that this example has limited alternatives of train headways (Table 3) and uses a conservative low bound. Comparing the current and optimized timetables, we find that the average proportions of energy consumption caused by passenger weights on all tracks are both around 23%. Although the proportions hardly change after optimization, there is a significant reduction in total energy consumption, implying that timetabling and the passenger assignment process affect energy consumption.
Overall, the solution algorithm produces satisfactory timetable solutions based on energy consumption reductions. The iterative process of EMSA is illustrated in Fig. 11, where EMSA-1, EMSA-2, EMSA-3 correspond to Eqs. (43), (44), and (45) for updating the impedances of running arcs, respectively. When the iterations get larger, the flow adjustment process is stable and convergent. The results show that Eqs. (44)-(45) (curves in red and green) accelerate the convergence. The comparisons between the original timetables and optimized timetable (when κ t = 1.3) are shown in Figs. 12a-b, in which Lines 1 and 2 are taken for demonstration. In the optimized timetable, all running trains of Line 1 incur a slight delay despite no change in the headway, whereas trains of Line 2 have both slightly increased headway and running times. The results are in line with the basic principles of energy consumption. Specifically, a larger headway means fewer trains on the planning horizon and a delay means a longer running time with lower energy consumption. Therefore, on the condition of transporting all passenger demand, the optimized timetable effectively reduces energy consumption.
To further show which OD pairs may change their paths under the optimized timetable, we record the OD pairs that have different least-cost paths from those under the original timetable. It is found that 896 OD pairs among the total 6588 OD pairs that need transfers have different least-cost paths. The visualization of origins and destinations of those pairs are shown in Fig. 13. The histogram height placed on the stations indicates the number of origins or destinations changing the least-cost path. We can see that the origins and destinations in the city center and northern part are more likely to change their paths under the optimized timetable. The maximum numbers of origins and destinations changing the least-cost path are 33 and 27, which are the Xingzhengzhongxin and Beidajie stations, respectively. This information is useful for aiding the operator to determine if the changes are intended.
A side effect is that passenger travel times may increase unavoidably. However, as shown in Fig. 14, compared with the original timetable, the largest increase of travel time is less than 18% of all OD pairs between the 94 stations of the URT network. It is found that the average travel time of all OD pairs is increased by 9.8% (the average travel time is 34.4 min with the original timetable). All setups being equal, energy-oriented timetabling cannot avoid the expenses of passenger travel time for reducing energy consumption (Chevrier et al., 2013;Yang et al., 2020). However, it appears acceptable to have considerably reduced energy consumption with a slight increase in passenger travel time in an initiative to combat environmental concerns (Mo et al., 2019a). Empirical research is required to achieve an appropriate trade-off between energy consumption and passenger travel time, which is however beyond the scope of this study.
Passenger volumes of Lines 1-4 are displayed in Fig. 15(a)-(d). We find that the passenger volume has two-peaks concentrating on certain tracks near the transfer stations. To transport all passengers with a periodic timetable, the peak volumes determine the frequency of the URT lines. Given that the train mass is much larger than the total passenger weight, a train itself accounts for a significant share of energy consumption. Therefore, a higher frequency means a larger fleet size and higher energy consumption. In Fig. 15(a)-(d), the obtained frequencies are denoted by the red dashed lines. To cut down energy consumption, a slight reduction in the peak volumes can effectively reduce the train frequencies to the green dashed lines, which correspond to the first larger alternative headways than the obtained ones. A larger value of κ t makes passengers less concentrate on the track near the transfer stations and results in a lower peak passenger volume. This also explains why a higher value of κ t leads to a smaller (better) upper bound.

Case 2.1: Effects of transferinterrupt of transfer opportunities
To show the effects of transfer between URT lines on the total energy consumption, we deliberately make an interrupt on transfer arcs at each transfer station (all transfer directions are interrupted at the station). Note that the interrupts on the transfer arcs at one station do not prevent the URT lines and the other transfer stations from functioning. The setup allows us to evaluate the influence of transfer opportunities due to potential emergencies on a single line separately. First, we consider the original passenger demand as in Case 1. When there is a separate interrupt at each transfer station, as shown in Fig. 16, it is found that the original passenger demand exceeds the capacity of the downgraded URT network and not all passengers can be transported. To show the influence after the interrupt, we define four indicators, "Energy consumption after an interrupt" (ECI), "Energy consumption increase after an interrupt" (EII), "Percentage of energy consumption increase after an interrupt" (EIP), "Maximum allowed percentage of the original passenger demand" (MAP). ECI, EII, and EIP refer to energy consumption change, while MAP reflects the vulnerability of the URT network. The ECI, EII, EIP, and MAP after the interrupt at each transfer station are shown in Fig. 17. It shows that the vulnerabilities of the transfer stations, indicated by MAP, are different. Amongst, the Xiaozhai is the weakest. If the transfer is not allowed, the network cannot even afford 60% of the original passenger demand. The interrupt at Wulukou has the least negative effect (− 5%) on the URT network. Since the Wulukou station is located in the central area, the interruption of transfer can be compensated by other transfer stations.
Furthermore, we use 50% of the original passenger demand (less than the values of MAP of all stations) proportionally for all ODs to show the increase in energy consumption. With the reduced demand, the optimal energy consumption is 36389.76 kWh before an interrupt at one of the transfer stations. Based on the indicator of ECI, we find that interrupts significantly increase energy consumption. Particularly, note that the EII and EIP may not change at the same pace with MAP. For example, the interrupts at the Xinzhenzhongxin and Wulukou stations have EII over 7000 kWh, which accounts for EIP about 20%. Interrupts at other stations have also caused about 10% in EIP. However, the Xinzhenzhongxin station has the least vulnerability for transfer, but it has the highest EIP. Combing the above results, from a reversed perspective, we find that transfer stations play an essential role in energy-saving.

Case 2.2: Effects of transferthe addition of transfer opportunities
To complement Case 2.1, we add transfer arcs based on Case 1 to further verify the influence of transfer opportunities on energy consumption. As shown in Fig. 18, two URT lines (Lines 5 and 6), which will come into operation soon 1 , have transfer opportunities with the existing lines at six stations. Since the detailed designs of the newly added URT lines are not available, we use the neighboring parallel sections to approximate the running times and energy consumption. To make interesting comparisons, we also assume that the added transfer opportunities (delineated within the two ellipses) would not introduce new passenger demand. Corresponding to Cases 1 and 2.1, we re-run the model under two scenarios and present the results below.
First, with the original passenger demand in Case 1, the energy consumption after separately adding transfer arcs of Lines 5 and 6 are 67945.27 kWh and 67613.02 kWh, respectively, meaning 2.3% and 2.8% of energy reductions, compared with the optimal solution (69537.97 kWh) in Case 1. When simultaneously adding transfer arcs of Lines 5 and 6, a higher energy reduction of 3.4% (67271.27 kWh) is achieved. Therefore, the transfer and synchronous optimization of multiple lines have a positive reduction in energy consumption. Second, considering 50% of the passenger demand in Case 1, we find that the vulnerability and energy consumption after an interrupt at the same station are reduced, compared to Case 2.1. In Fig. 19, '-5 ′ , '-6 ′ , and '-5&6 ′ are attached to the indicators referring to the addition of transfer arcs with Line 5 only, Line 6 only, and Lines 5-6, respectively. As seen in the histogram, the EII values after the interrupt of transfer opportunities associated with the Beidaijie, Wulikou, Xiaozhai, and Dayanta stations are negative. The EIP values after adding the transfer opportunities show obvious reductions with reference to the green dashed curve (before the addition). Moreover, the vulnerability of the URT network has been reduced after the addition of transfer opportunities, as shown in Fig. 20. The original MAP, denoted by the purple dashed curve, is lower than the solid curves. After adding transfer arcs of Line 5 or Line 6, the MAP values are different due to the interrupt at different stations. Even when the transfer arcs at Beidajie or Wulukou are interrupted, the network can still afford all (100%) of the original passenger demand. When the interrupt happens at Tonghuamen, Xiaozhai, or Danyanta, the maximum affordability of the network is significantly improved after adding the transfer opportunities. Only for the interrupt at the Xinzhengzhongxin station (located in the suburban area), the added transfer opportunities have no effect on the MAP.
Overall, through the above numerical analyses, it can be concluded that the proposed model and the solution algorithm perform well for generating energy-efficient timetables of the URT network with multiple lines. The results also show that the transfer opportunities and passenger path choice significantly influence the total energy consumption of the URT network. For the URT operator, (EII: Energy consumption increase after an interrupt; EIP: Percentage of energy consumption increase after an interrupt)  the transfer design and control are not only crucial to vulnerability reduction but also important for the energy-saving of the URT network.

Conclusions and future work
Existing energy-efficient timetabling has been limited to single URT lines. In the paper, we proposed a bi-level model framework that formulates energy-efficient timetabling of an entire URT network as MILP in the UL and UE-based path choice responding to a timetable in the LL. A heuristics algorithm involving a relaxation process and an MSA module was developed to find near-optimal timetable solutions. The efficiency and effectiveness of the proposed model were demonstrated in experimental examples considering the URT network of Xi'an (China) as a study area. The proposed model can achieve a significant reduction (26.4%) of energy consumption with the compromise of a minor (9.8%) increase in the average travel time compared to the current timetable. To highlight the influence of transfer opportunities and timetable synchronization, sensitivity analyses based on the deletion and addition of transfer arcs are conducted. The results show that transfer opportunities and synchronization between multiple URT lines contribute to saving energy and improving the resilience of the URT network. These results confirm the inaccuracy and insufficiency of the existing energy-efficient timetabling models that ignore the interconnections between multiple URT lines.
Moreover, this study paves the way to several promising research directions. First, the trade-off between energy consumption and passenger travel time should be considered, for which multi-objective energy-efficient timetabling is worth investigation. Second, the solution algorithm should consider space-time networks in the full temporal dimension. In theory, higher resolutions improve the fidelity of the energy-efficient tabling but require more computation resources. Therefore, speedup techniques should be developed for the solution algorithm. Third, some new URT technologies, such as energy regenerative and energy storage devices, could be incorporated into the model framework. Fourth, more realistic passenger path choice behavior should be incorporated in the URT network. Finally, it is important to investigate the relationship between the topology of the URT networks and energy consumption for designing the URT networks. We will address these issues in our future work.

Table A1
The running time and energy consumption of speed profile in three levels.