Hierarchical framework for coordinate optimization of train timetable and electric motor train unit circulation plan

Abstract

To meet the increasing demand for improving railway service quality while using the precious mobile resources reasonably and economically, this research proposes a hierarchical approach that integrates the model for constructing an efficient electric motor train unit circulation plan into the model for optimized timetable design without predefining timetable details. A simulated annealing algorithm for solving the timetabling (main) model is designed, in which the neighborhood system is pretreated. A special tree construction–based branch-and-bound algorithm is improved for solving the electric motor train unit circulation planning (sub)model. The results of the numerical experiment verify the effectiveness of the proposed method. Specifically, compared to the randomly generated initial solution, the optimal solution obtained by the proposed methods reduces the total travel time by 686 min and reduces the number of electric motor train units by 5. The number of electric motor train units needed by the proposed method is on average at least two less than the method that handles the problem in a sequential way. Railway operators can implement this approach for balancing the efficiency of timetable and the quality of electric motor train unit circulation plan within a reasonable computation time when scheduling trains in railways.

Keywords

High-speed railway network train timetabling electric motor train unit circulation plan coordination optimization

Introduction

Train timetable and electric motor train unit (EMU) circulation plan are among the most important plans in the operation of high-speed railway systems. In the timetabling problem, the exact departure and arrival times of each train at each station along its route are determined. Certain operational rules such as track capacity limitation and achievement of the highest possible passenger service quality (e.g. travel time) should be considered. A number of studies related to this problem are available in various literatures. Caprara et al.¹ and Brännlund et al.² developed the Lagrangian relaxation–based heuristics to compute a timetable. The track capacity constraints were relaxed in order to make the problem to be handled more easily. Zhou and Zhong³ proposed two solving algorithms: a branch-and-bound algorithm and a beam search algorithm to deal with the double-track train-scheduling problem with multiple objectives. Lee and Chen⁴ presented an optimization heuristics to solve the train timetabling problem when considering the problems of train pathing. The method was able to deal with real-sized instances. Castillo et al.⁵ proposed an optimization method to solve the train-scheduling problem for a double- and single-tracked bi-directional railway network. Some excellent techniques were incorporated into the method to speed up the optimization solution of the problem. Cacchiani et al.⁶ modeled the problem of timetabling at a station. The model did well in considering both complex rules for correct junction traversal, and constraints on the maximum number of trains allowed for the simultaneous presence at a station. Chu et al.⁷ studied the bi-objective optimal urban rail train scheduling. The presented method aimed at finding a good train timetable which achieves considerable service quality with as little operation cost as possible.

The EMU circulation planning problem needs to determine how to allocate the EMUs to the train trips with a given timetable under the condition that the EMUs should be maintained properly. Specifically, the circulation plan consists of a set of connection relationships between train trips. These connection relationships indicate which train trips each EMU should fulfill in each day and when/where to arrange the first-level maintenance for the EMUs. All of the EMUs are operated in a periodic manner according to the circulation plan. In railways, there are often several interrelated train routes. Each route corresponds to a part of railway line. It is more reasonable to utilize the EMUs throughout the interrelated routes than separating the EMUs into several groups and each for a different route (see, for example, Zhao et al.⁸).

Due to the complexity of the timetabling problem itself, the EMU circulation plan is often ignored when finding an optimal timetable for trains. The computations of these two plans are handled in a sequential way. We refer to this method as the ordinary method in this article. The timetable is served as the input when finding EMU circulation plan. Peeters and Kroon⁹ described a branch-and-price algorithm to deal with the problem of efficient circulation of train units within a certain scope of railway line given the timetable and the passengers’ seat demand. Alfieri et al.¹⁰ presented a solution approach to determine the appropriate numbers of train units of different types together with their efficient circulation on a single line. The approach is based on an integer multi-commodity flow model. Cacchiani et al.¹¹ proposed a heuristics for the train unit assignment problem. The proposed method is based on the natural Lagrangian relaxation of a natural integer linear programming and turns out to be much faster in practice and still providing solutions of good quality. Abbink et al.¹² considered the train types and subtypes when dealing with the problem of allocating units of rolling stock to the train series, such that as many people as possible can be transported with a seat, especially during the rush hours. Shi et al.¹³ designed a simulated annealing algorithm for the EMU circulation planning problem. The algorithm is based on the penalty function and 3-opt neighborhood structure for the circular permutation of all trains. Zhao and Tomii¹⁴ introduced a probability-based local search algorithm for solving the connection of the trains and the generation of the maintenances of the problem after it was transformed into the traveling salesperson problem.

Besides the deterministic models for strategic scheduling of trains or circulation planning of the EMUs, there are also methods that focus on designing timetables or EMU circulation plans to handle the different types of disturbances, such as delays or any other unforeseen events occurring at an operational level.^15,16 These methods are mainly based on robustness optimization or delay (disturbance) control. YL Wei et al.¹⁷ proposed a powerful reliable controller to preserve system stability, accommodate the component faults, and maintain the preferable system performances using Markov chain to describe component failures commonly occurred in dynamical (fuzzy) systems. The various controlling approaches based on Markov models are able to handle different types of dynamic systems subjected to random abrupt changes in their structures.^18–21

Since the timetable provides the premises for the EMU circulation plan, even a subtle change in the departure and arrival times of some of the trains would have a direct influence on the potential optimal EMU circulation plan. When the EMU circulation plan is found, the number of EMUs needed is determined accordingly. As shown in Figure 1, although the numbers of train trips in the timetable (a) and (b) are the same, the structures of these two timetables are different. As a result, the circulation plan in (a) consists of at least four circulation fragments (①-④). The circulation plan in (b) consists of at least five circulation fragments (①-⑤). A circulation fragment represents the train trips, connections of train trips, and first-level maintenance that an EMU should fulfill in one period (e.g. 24 h) of schedule. The whole circulation plan is formed by connecting all the fragments through the overnight connections. The number of the necessary EMUs is equal to the number of circulation fragments contained in a circulation plan.

Figure 1.

Circulation fragments in two different circulation plans.

Therefore, the number of EMUs needed for fulfilling all the tasks of the train trips is not only related to how to build connection relationships between these train trips but also related to what the exact departure and arrival times of each train trip have been. It is reasonable to consider the EMU circulation plan when calculating the timetable for the train trips. The existing studies did not put enough emphasis on making a coordination optimization of the detailed timetable and EMU circulation plan. The particular interest of this article is to find an optimal timetable which is beneficial to efficiently utilize the EMUs. The main contribution of this article is that we propose an approach to optimize the timetable and EMU circulation plan in a unified way. The proposed hierarchical approach integrates the construction of an efficient EMU circulation plan into the optimized design of the timetable under the condition that the timetable details (all departure and arrival times, train order, the stations for overtaking operations, etc.) are not predefined. Furthermore, the algorithms for solving the main model and the sub model are efficiently designed.

The structure of this article is organized as follows: section “Problem illustrations” first describes the conceptual overview of the hierarchical framework proposed for the coordinate optimization of timetable and EMU circulation plan. Then, the space-time-network-based problem illustration is presented. Section “Formulation of the coordinate optimization framework” first defines the notations of the parameters and variables used in the proposed framework. Then, the main timetabling model and the EMU circulation planning submodel are formulated. Section “Solution method for the coordinate optimization framework” describes the procedure of the algorithm for solving the main model which is incorporated with the algorithm for solving the sub model. Section “Numerical example” presents the comparison of the results of the experiments with the proposed method and the corresponding methods in contrast. Section “Conclusion” summarizes the proposed method and mentions the future tasks of this research.

Problem illustrations

Conceptual hierarchical framework

There are two types of objectives in the problem of coordinate optimization of timetable and EMU circulation plan. First, a timetable with optimal efficiency should be found. Second, the quality of the EMU circulation plan calculated based on the found timetable is also optimal among those calculated based on other timetables. Although these two objectives are integrated in the single optimization problem, the variables used for the evaluation of the EMU circulation plan are separated from those variables used for the construction of the optimal timetable (as shown in Figure 2). The determination of these two types of variables is relatively independent. At the same time, the decision space of finding the optimal EMU circulation plan is determined by the decision space of finding the optimal timetable. Therefore, the computations in the proposed framework for coordinate optimization of the timetable and EMU circulation plan naturally consist of a timetabling model and an EMU circulation plan calculation model. The former is the main model and the latter is the submodel.

Figure 2.

Framework for coordinate optimization of the timetable and EMU circulation plan.

Specifically, the computation process of the framework could be summarized as follows: at each iteration of the main model, a feasible timetable is generated without predefined details; the quality of the circulation plan based on the iterative timetable is evaluated through the submodel; the total objective value of the iterative timetable is compared with the best timetable found so far, and the best timetable is possibly updated accordingly; Finally, the balanced optimal timetable and the EMU circulation plan are found.

Space-time-network-based problem description

For each train with a certain physical route on the railway network, a space–time network is associated with it. The space–time network of each train consists of several vertex sets and three types of directed arc sets. The position of each vertex set on the vertical axis corresponds to the station along its route. Each vertex set consists of linearly arranged vertices which are parallel to the horizontal axis and represent the discretized timestamps that the train could depart from or arrive at the station. The first type of arc is the virtual arc. All of the departure vertices at the first station and the arrival vertices at the last station are connected with an artificial source vertex $σ$ and an artificial terminal vertex $τ$ by the virtual arcs. The second type of arc is the section arc, which represents a possible running activity of the train from a departure vertex of a station to the arrival vertex of the next station. The third type of arc is the station arc, which represents dwelling activity of the train. Any arrival vertex at a station is connected with all possible departure vertices at the same station by the station arcs if the train is scheduled to dwell at that station. The directions of section arcs or station arcs are similar to the running direction of the train. Each arc is associated with a cost, which is equal to the duration between its start vertex and end vertex. The cost of the virtual arc is set to be zero. The timetable of a train can then be obtained by selecting a path from $σ$ to $τ$ in its space–time network.

Given a topology of a railway network where the railway is double track, a multi-graph could be constructed by combining the space–time network associated with each train to be scheduled (as shown in Figure 3). Let the connection arc denote the possible connection relationship between trains at a station. There are also train connection arcs in the multi-graph between each pairs of departure trains and arrival trains at stations such that the duration between the arrival time of a train and the departure time of another train is at least equal to the minimum connection time. The direction of a connection arc is from the arrival train to the departure train which means that after fulfilling the train trip and arriving at the station, the EMU will continue to fulfill the next train trip departing from the same station. The minimum connection duration of EMU at a station is the necessary criterion of the EMU dwelling time for performing the passengers’ alighting/boarding, EMU cleaning, and servicing activities before starting the next train trip.

Figure 3.

Space–time multi-graph of trains.

For clarity purposes, Figure 3 only shows the example of the selected path and some possible range of candidate paths associated with the selected path for each train. Then, the problem of the coordinate optimization of the timetable and the EMU circulation plan aims to find a feasible path for each train and a set of connection arcs linking all the trains in the space–time multi-graph under the condition that the efficiency of the timetable and the quality of the EMU circulation plan are balanced. In order to generate a feasible circulation plan, when selecting the connection arc, the cumulative running time/distance of EMU should be considered. The cumulative running time/distance of EMU could not exceed the upper bound. Otherwise, maintenance must be arranged (in advance) for the EMU at the stations linked with an EMU depot. In case of arranging maintenance only, the connection arcs whose durations are longer than the necessary maintenance time can be selected. For example, in Figure 3, stations 1 and 2 are linked with an EMU depot. The maintenance can be arranged at these stations in the circulation plan.

Formulation of the coordinate optimization framework

Notations of parameters and variables

The parameters include the following: K is the set of trains; $K'$ is the set of high-speed trains; $K''$ is the set of medium-speed trains; k or $k'$ is the index of train; $S_{k}$ is the set of stations along the route of train k; ${\bar{S}}_{k} = S_{k} - {f_{k}, t_{k}}$ , where $f_{k}, t_{k}$ is the first and last stations in $S_{k}$ ; $S'$ is the set of first stations of all the routes; s is the index of station; $g_{s, k}$ is the minimum time for train k to dwell at station s; $E_{k}$ is the set of arcs in the space–time network of train k; $E'_{k}$ is the station arcs in $E_{k}$ ; e, $e'$ , or $e''$ is the index of arc; $V_{k}$ is the set of vertices in the space–time network of train k; $N_{k, s}$ is the set of departure vertices at station s; $R_{k, s}$ is the set of arrival vertices at station s; v is the index of vertex; $c_{e}$ is the cost of the arc e; $T_{s}$ is the set of trains whose first station is station s; $δ_{s}$ is the mean value of the intervals of departure times when the interval of each adjacent train departure times at station s is the same during the planning horizon; $Θ_{v}^{+}$ is the set of arcs leaving vertex v; $Θ_{v}^{-}$ is the set of arcs entering vertex v; $B_{e}$ is the set of trains whose route goes through the section identified by the arc e; $U_{k, e}$ is the set of arcs in $E_{k}$ that are incompatible with arc e in the section; $d_{e}$ is the timestamp of the starting vertex of the arc e; $a_{e}$ is the timestamp of the ending vertex of the arc e; $η_{k}$ is the upper bound of the extra waiting time of medium-speed train k generated due to being overtaken by other high-speed trains; $π_{k}$ is the mileage from the first station to the last station of the route of train k; the binary parameter $φ_{k, k'}$ indicates whether there is a candidate connection relationship between train k and $k'$ ; $m_{e}$ and $h_{e}$ are the arrival train and departure train in the connection arc e; $α$ and $β$ are the upper bound of the allowable accumulated traveling mileage and traveling time for a EMU after a maintenance; $ε$ is the minimum duration required for the connection between trains; $\bar{ε}$ is the minimum duration required for the maintenance; and H is the duration of the planning horizon.

The decision variables include the following: the binary variable $x_{k, e}$ indicates whether the arc e is selected in the space–time path of train k; the binary variable $y_{k, k'}$ indicates whether the connection relationship is built between train k and $k'$ ; the binary variable $z_{k, k'}$ indicates whether the maintenance is arranged in the duration of the connection between train k and $k'$ ; The intermediate binary variable $ρ_{k, k'}$ indicates whether the duration of the connection between train k and $k'$ enables the arrangement of the maintenance; $l_{k, k'}$ indicates the duration of the connection between train k and $k'$ ; $Δ_{s, k}$ is the interval between the departure time of train k (except for the last train) and the departure time of the train behind train k at station s; intermediate variable $P_{k}$ is the set of trains from each of which there is a path to the train k in the multi-graph when arbitrarily deleting two connection arcs in the whole circulation plan.

Formulation of the main model

The main model in the framework selects the set of space–time paths for each train under the condition of satisfying the safety rules. The timetable constructed based on the selected paths should have an optimal efficiency. The optimal quality of the circulation plan is another objective used to evaluate the obtained timetable. But the detailed selections of the train connection arcs and calculation of the objective value of the optimal circulation plan are realized through the solving of the submodel. Specifically, the main model is formulated as follows

\begin{matrix} min μ_{1} \sum_{k \in K} \sum_{e \in E_{k}} c_{e} x_{k, e} \\ + μ_{2} \sum_{s \in S'} \sum_{k \in T_{s}} {(Δ_{s, k} - δ_{s})}^{2} + μ_{3} C (L, Y, Z) \end{matrix}

(1)

Subjected to

\sum_{e \in Θ_{σ_{k}}^{+}} x_{k, e} = 1 \forall k \in K

(2)

\sum_{e \in Θ_{v}^{+}} x_{k, e} = \sum_{e \in Θ_{v}^{-}} x_{k, e} \forall k \in K, v \in V_{k}, s \in S_{k}

(3)

\sum_{k \in B_{e'}} \sum_{e \in U_{k, e'}} x_{k, e} \leq 1 \forall e' \in E

(4)

\sum_{v \in N_{k, s}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k, e} - \sum_{v \in R_{k, s}} \sum_{e \in Θ_{v}^{-}} a_{e} x_{k, e} \geq g_{s, k} \forall s \in {\bar{S}}_{k}, k \in K

(5)

\sum_{e \in E_{k}^{''}} c_{e} x_{k, e} \leq η_{k} \forall k \in K''

(6)

Objective function (1) contains three components. The first component is the total travel time. The reduction in travel time would improve the service quality of the railway enterprise and the satisfaction of passengers. The second component is the total mean variance of intervals of each adjacent train departure times at each station $s \in S'$ . By incorporating this component, we hope to obtain a set of departure times such that each departure time interval deviates from the mean value as little as possible. The balancing of the departure time intervals of trains is beneficial to use the fixed facilities of stations more efficiently and to improve the stability of the timetable. The third component $C (L, Y, Z)$ is the cost of the EMU circulation plan. A high quality of EMU circulation plan with minimum cost would save the number of EMUs needed to fulfill all the train trips. $C (L, Y, Z)$ will be computed through the solving of the submodel. L, Y, and Z are the vectors of the variables $l_{k, k'}$ , $y_{k, k'}$ , and $z_{k, k'}$ , respectively. $l_{k, k'}$ is an intermediate variable associated with $x_{k, e}$ : if $\sum_{e \in Θ_{σ_{k'}}^{+}} a_{e} x_{k', e} - \sum_{e \in Θ_{τ_{k}}^{-}} d_{e} x_{k, e} \geq ε$ , then

l_{k, k'} = \sum_{e \in Θ_{σ_{k'}}^{+}} a_{e} x_{k', e} - \sum_{e \in Θ_{τ_{k}}^{-}} d_{e} x_{k, e}

(7)

Else

l_{k, k'} = \sum_{e \in Θ_{σ_{k'}}^{+}} a_{e} x_{k', e} - \sum_{e \in Θ_{τ_{k}}^{-}} d_{e} x_{k, e} + H

(8)

Constraint (2) indicates that only one arc can be selected among those leaving the virtual starting vertex for each train. Constraint (3) is the flow balance constraints, which ensures to select a path for each train in its space–time network. Constraint (4) imposes that there are no conflicts between the space–time paths of trains. Constraint (5) imposes that the dwelling time of each train at each station should not be less than the lower bound. Constraint (6) imposes that the total extra waiting time of medium-speed trains due to being overtaken by high-speed trains should be within a certain scope. The extra waiting time is related to the utilization of the capacity of the route and can be determined through experimentation. This group of constraints could also be used as a supplement to the first component of the objective.

Formulation of the submodel

On the basis of the variables $l_{k, k'}$ that are transmitted from the main model, the submodel mainly makes the decisions of constructing the connection relationships between trains and arranging the maintenance tasks in the whole circulation plan. The purpose of the submodel is to find an optimal circulation plan when given each set of the variables $l_{k, k'}$ . The submodel could be formulated as follows

C (L, Y, Z) = min μ_{4} \sum_{k \in K} \sum_{k' \in K} l_{k, k'} y_{k, k'} + μ_{5} \sum_{k \in K} \sum_{k' \in K} z_{k, k'}

(9)

Subjected to

\sum_{k \in K} y_{k, k'} = 1 \forall k' \in K

(10)

\sum_{k' \in K} y_{k, k'} = 1 \forall k \in K

(11)

y_{k, k'} \leq φ_{k, k'} \forall k, k' \in K

(12)

\sum_{e \in F_{e', e''}} l_{m_{e}, h_{e}} y_{m_{e}, h_{e}} + \sum_{e \in F_{e', e''}} θ_{h_{e}} \leq α \forall e', e'' \in W

(13)

\sum_{e \in F_{e', e''}} π_{h_{e}} \leq β \forall e', e'' \in W

(14)

z_{k, k'} \leq ρ_{k, k'} y_{k, k'} \forall k, k^{'} \in K

(15)

P_{k} + P_{k'} = | K | \forall k, k^{'} \in K

(16)

\sum_{k \in K} \sum_{k' \in K} z_{k, k'} \geq 1

(17)

The objective function (9) includes two components. The first component is the total connection time of all the trains. Then, the total number of EMUs $ρ (X, Y)$ needed to fulfill the tasks of all the train trips with given variables $x_{k, e}$ and $y_{k, k'}$ could be computed as follows

ρ (X, Y) = \frac{(\sum_{k \in K} \sum_{e \in E_{k}} c_{e} x_{k, e} + \sum_{k \in K} \sum_{k' \in K} l_{k, k'} y_{k, k'})}{H}

(18)

The second component is the times of the maintenance in the circulation plan. Constraints (10) and (11) impose that there is only one connection arc entering and leaving each train. Constraint (12) imposes that the connection arcs should be selected from the candidate arcs between trains. Constraints (13) and (14) impose that the accumulated traveling mileage and traveling time for an EMU after maintenance should be less than the corresponding upper bound. Intermediate variable W indicates the set of maintenance arcs. Intermediate variable $F_{e, e'}$ indicates the set of connection arcs between each adjacent maintenance arcs $e, e'$ in the whole circulation plan. Intermediate variable $θ_{k}$ indicates the travel time of train k. Constraint (15) imposes that the connection arc could be arranged for maintenance only if the arc satisfies the corresponding maintenance condition. Constraint (16) avoids the subloop in the circulation plan. Constraint (17) ensures that there is at least one maintenance arc in the circulation plan.

Solution method for the coordinate optimization framework

It can be seen from each of the model that they are both complicated integer programming problem. In order to ensure the speed in solving the proposed framework, which integrates these two models, a simulated annealing algorithm is designed. It is embedded with an improved branch-and-bound method. The neighborhood of the tentative solution when solving the main model in the simulated annealing process is predetermined, although the departure time of each train could be within the whole planning horizon. In fact, the experiment will show that this neighborhood designing mechanism has a great impact on ensuring the quality of the obtained solution. The procedure of the algorithm is as follows:

Step 1: initialization—set the initial temperature $T > 0$ , maximum times of internal iteration $j_{\max} > 0$ , and the current internal iteration $j = 0$ . Let $\bar{ξ}$ denote a complete random solution of the framework. Set the initial optimal value of the objective function $F (ξ_{min}) = \infty$ .

Step 2: according to the number of elements in the set of high-speed trains $T'_{s}$ and the set of medium-speed trains $T {''}_{s}$ , where $T_{s} = T {''}_{s} \cup T'_{s}$ , generate a random vector $(r_{1}^{j}, r_{2}^{j}, \dots, r_{| {\bar{T}}_{s} |}^{j}), {\bar{T}}_{s} \in {T {''}_{s}, T'_{s}}$ for each of them. The value of the element $r_{k}^{j}$ is as follows

r_{k}^{j} = \frac{sign (u_{k}) (χ_{k} - ϖ_{k})}{| {\bar{T}}_{s} |}, k \in {\bar{T}}_{s}

(19)

where $u_{1}, u_{2}, \dots, u_{k}$ are a group of stochastic variables that are uniformly distributed in the interval $[- 1, 1]$ and independent from each other. $sign (\cdot)$ is the sign function. $χ_{k}$ and $ϖ_{k}$ are the earliest and the latest departure time in H allowed for train k.

Step 2.1: for the tth train $k \in T'_{s}$ , select the virtual arc entering the first station as follows

\sum_{e \in Θ_{σ_{k}}^{+}} a_{e} x_{k, e} = r_{k}^{j} + \frac{t (χ_{k} - ϖ_{k})}{| {T'}_{s} |}

(20)

Since the high-speed train has the highest rank and could not be influenced by other trains, the remaining part of the path of train $k \in T'_{s}$ in the space–time network is determined using the shortest path method.

Step 2.2: for each train $k \in T {''}_{s}$ , set the index of the current section ${\hat{s}}_{k} = 1$ . The departure vertex in the first station is selected the same as equation (20) but the remaining part of the path of train $k \in T {''}_{s}$ is determined according to the following path-finding procedure: set the temporary variable $ω_{1} = \sum_{e \in Θ_{σ_{k}}^{+}} a_{e} x_{k, e} + γ_{k, 1} + t_{k, q} + t_{k, t}$ , where $γ_{k, 1}$ indicates the pure running time of train k in the first section. $t_{k, q}$ and $t_{k, t}$ are the additional time for departure and arrival of train k. Update the current section ${\hat{s}}_{k} = {\hat{s}}_{k} + 1$ .

Step 2.3: if the current section ${\hat{s}}_{k}$ is not the last section in the route of train k, set the temporary variable ${\bar{ω}}_{2} = ω_{2} = ω_{1} + g_{s', k}$ , where $s'$ is the entering station of section ${\hat{s}}_{k}$ . Set $k'$ to be the first high-speed train in section ${\hat{s}}_{k}$ . Go to step 2.4. Otherwise, go to step 2.8.

Step 2.4: if $k'$ is not the last high-speed train in the section ${\hat{s}}_{k}$ and if one of the following three conditions is satisfied (as shown in Figure 4), go to step 2.5, otherwise, set $k'$ to be the next high-speed train in section ${\hat{s}}_{k}$ and return to step 2.4. If $k'$ is already the last high-speed train in the section ${\hat{s}}_{k}$ , go to step 2.6.

(a) The interval between $ω_{2}$ and the departure time of train $k'$ at station $s'$ is less than the minimum safety headway I

| \sum_{v \in N_{k', s'}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k', e} - ω_{2} | < I

(21)

(b) The interval between the arrival time of train k and the arrival time of train $k'$ at station $s''$ is less than the minimum safety headway I, where $s''$ is the leaving station of section ${\hat{s}}_{k}$

| ω_{2} + γ_{k, {\hat{s}}_{k}} + t_{k, q} + t_{k, t} - \sum_{v \in R_{k', s''}} \sum_{e \in Θ_{v}^{-}} a_{e} x_{k', e} | < I

(22)

Figure 4.

Three conditions for the overtaking.

In condition (b), since train k is currently not sure whether it needs to dwell at station $s''$ for the reason of being scheduled or overtaken by the high-speed trains including $k'$ , the arrival safety headway checking is based on the condition that train k will dwell at station $s''$ , that is, $t_{k, q}$ and $t_{k, t}$ should be considered.

\begin{matrix} (ω_{2 -} \sum_{v \in N_{k', s'}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k', e}) \\ (ω_{2} + γ_{k, {\hat{s}}_{k}} + t_{k, q} + t_{k, t} - \sum_{v \in R_{k', s''}} \sum_{e \in Θ_{v}^{-}} a_{e} x_{k', e}) < 0 \end{matrix}

(23)

Step 2.5: update the temporary variable $ω_{2} = \sum_{v \in N_{k', s'}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k', e} + I$ . It means that train k is overtaken by train $k'$ at station $s'$ . Update the $k'$ to be the next high-speed train in section ${\hat{s}}_{k}$ . Go to step 2.4.

Step 2.6: if ${\bar{ω}}_{2} = ω_{2}$ , the actual arrival vertex of the arc to be selected for train k in section ${\hat{s}}_{k} - 1$ is updated as

\sum_{v \in R_{k, s - 1^{''}}} \sum_{e \in Θ_{v}^{-}} a_{e} x_{k, e} = ω_{1} - (1 - ψ_{k, s - 1^{''}}) t_{k, t}

(24)

where $s - 1''$ is the leaving station of section ${\hat{s}}_{k} - 1$ , $ψ_{k, s - 1''}$ is the binary parameter indicating whether train k is planned to dwell at station $s - 1''$ . Besides, the actual departure vertex of the arc to be selected for train k in section ${\hat{s}}_{k}$ is updated as

\sum_{v \in N_{k, s'}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k, e} = \sum_{v \in R_{k, s - 1^{''}}} \sum_{e \in Θ_{v}^{-}} a_{e} x_{k, e} + g_{s', k}

(25)

Otherwise, the corresponding arrival vertex and departure vertex of the arc to be selected for train k are determined as follows, and go to step 2.7

\sum_{v \in R_{k, s - 1^{''}}} \sum_{e \in Θ_{v}^{-}} a_{e} x_{k, e} = ω_{1}

(26)

\sum_{v \in N_{k, s'}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k, e} = ω_{2}

(27)

Step 2.7: update the temporary variable $ω_{1} = \sum_{v \in N_{k, s'}} \sum_{e \in Θ_{v}^{+}} d_{e} x_{k, e} + γ_{k, {\hat{s}}_{k}} + t_{k, q} + t_{k, t}$ and update ${\hat{s}}_{k} = {\hat{s}}_{k} + 1$ . Go to step 2.3.

Step 2.8: if the whole space–time path of train $k \in T {''}_{s}$ has no conflicts with other trains in $T {''}_{s}$ whose space–time paths in the tentative point of the algorithm have already been determined, go to step 3. Otherwise, update the virtual arc entering the first station as follows

\sum_{e \in Θ_{σ_{k}}^{+}} a_{e} x_{k, e} = r_{k}^{j} + \frac{t (χ_{k} - ϖ_{k})}{| {T'}_{s} |} + (- 1)^{ε_{j}} ε_{j}

(28)

where $ϵ_{j} = 1$ , set $ϵ_{j} = ϵ_{j} + 1$ , and go to step 2.2 (reconstruct the space–time path for train k), until train k has no conflicts with other trains in $T {''}_{s}$ , go to step 3.

If the space–time networks of trains in the multi-graph are abstracted into another kind of vertices, then the multi-graph could be transformed into a connection network. Figure 5 shows the connection network transformed from the multi-graph shown in Figure 3. Given a network with $n + 1$ vertices, a 1-tree is defined as a set of $n + 1$ arcs that span the graph (see, for example, Fisher.²²). The Degree Constraint Minimum 1-Tree (DCM1-T) is a 1-tree such that the total weights of the set of the arcs in the 1-tree are minimum and that the degree of each node does not exceed the limitation of the degree constraint.

Figure 5.

Connection network transformed from the space–time multi-graph.

It can be seen from Figure 5 that the problem of calculating the EMU circulation plan is equal to the problem of finding a DCM1-T in the connection network such that the in-degree and out-degree of each vertex in the DCM1-T are equal to 1 and that the maintenance is taken into consideration.

Step 3: construct the connection network of the set of the train K. Set the maximum times of iteration of calculating the EMU circulation plan to be $d_{max}$ , and the current iteration $d = 0$ . Set the initial upper bound of the EMU circulation plan $Φ = \infty$ . Set the initial branching level $i = 0$ , the initial index of branch (at level i) $m (i) = 0$ , the initial right lower bound $R (i) = \infty$ , and the initial number of branches to be $λ$ . We define $E N_{m (i)}$ and $E X_{m (i)}$ as the set of connection arcs included and excluded in the circulation plan calculated by the branch m of branching level i. Obviously, these two sets are used to fix variables in the branching process. Set the initial set $E N_{0 (1)}$ to include an arbitrary connection arc which could be used as the maintenance arc and $E X_{0 (1)} = \emptyset$ .

Step 3.1: if $d < d_{max}$ , call the Degree Constraint Minimum 1-Tree Generating (DCM1-TG) algorithm²³ on the basis of the set $E N_{m (i)}$ and the circulation plan $C_{m (i)}$ is obtained. Set $i = i + 1$ and go to step 3.2. Otherwise, go to step 4.

Since the problem of finding DCM1-T is NP-complete (see, for example, Alexandre and Abilio²⁴), the DCM1-TG algorithm is able to find an approximate lower bound for the EMU circulation plan. The circulation plan obtained by the DCM1-TG algorithm could not be feasible because the maintenance is not considered. In order to find an optimal feasible circulation plan, an improved branch-and-bound method is applied as follows

Step 3.2: according to the set of $E N_{m (i - 1)}$ and $E X_{m (i - 1)}$ , determine the set $E N_{λ (i)}$ and $E X_{λ (i)}$ as follows

E N_{λ (i)} = E N_{m (i - 1)} \cup ℓ \overset{{{\bar{e}}_{λ (i)}}}{\cup}

(29)

E X_{λ (i)} = E X_{m (i - 1)} \cup^{} {{\bar{e}}_{1 (i)}, {\bar{e}}_{2 (i)}, \dots, {\bar{e}}_{λ - 1 (i)}}

(30)

where $ℓ$ is the maximum fragment in circulation plan $C_{m (i)}$ that starts from the last arc (sorted by the direction of the arcs) in $E N_{m (i - 1)}$ and satisfies that the accumulative running distance and the running time do not exceed the allowable upper bounds. ${{\bar{e}}_{1 (i)}, {\bar{e}}_{2 (i)}, \dots, {\bar{e}}_{λ (i)}}$ is a set of connection arcs, each of which starts from the arrival train of the last arc in $ℓ$ and is able to be used as the maintenance arc.

Call the DCM1-TG algorithm on the basis of the set $E N_{λ (i)}$ and $E X_{λ (i)}$ and the circulation plan $C_{λ (i)}$ is obtained. Set the left lower bound $L (i) = F' (C_{λ (i)})$ , where $F' (C_{λ (i)})$ is the objective value of the submodel calculated with the circulation plan $C_{λ (i)}$ . Update $m (i) = λ (i) - 1$ . Determine the set $E N_{m (i)}$ and $E X_{m (i)}$ . Call the DCM1-TG algorithm on the basis of the set $E N_{m (i)}$ and $E X_{m (i)}$ and the circulation plan $C_{m (i)}$ is obtained. Update the right lower bound $R (i) = F' (C_{m (i)})$ .

Step 3.3: if $L (i) < Φ$ , go to step 3.4, otherwise, go to step 3.5.

Step 3.4: if the circulation plan corresponding to $L (i)$ is feasible with respect to the maintenance constraints, update $Φ = L (i)$ and go to step 3.5, otherwise, update $i = i + 1$ and go to step 3.2.

Step 3.5: if $R (i) < Φ$ , update left lower bound $L (i) = R (i)$ and go to step 3.6, otherwise, go directly to step 3.6.

Step 3.6: Update $m (i) = m (i) - 1$ . If $m (i) = 0$ , go to step 3.7, otherwise, if $m (i) > 1$ , determine the set $E N_{m (i)}$ and $E X_{m (i)}$ . Call the DCM1-TG algorithm based on the set $E N_{m (i)}$ and $E X_{m (i)}$ and the circulation plan $C_{m (i)}$ is obtained. Update $R (i) = F' (C_{m (i)})$ . If $m (i) = 1$ , update $R (i) = \infty$ . After the $R (i)$ is determined, check whether $R (i) > Φ$ and go back to step 3.6, otherwise, go to step 3.3.

Step 3.7: backtrack to the previous level—update $i = i - 1$ ; if $i > 1$ , go to step 3.5, otherwise update $λ = λ + 1$ . If $λ < λ_{max}$ , go to step 3.1. Otherwise, go to step 3.8.

Step 3.8: update $d = d + 1$ , $i = 0$ . Randomly select a feasible part of circulation plan generated at the iteration $d - 1$ and replace the last maintenance arc of the selected part of the circulation plan, or randomly select a maintenance arc in the circulation plan generated at the iteration $d - 1$ . Initialize the $E N_{0 (1)}$ as the set of arcs included in the selected part of the circulation plan or as the selected maintenance arc. Go to step 3.1.

Step 4: a new complete random solution of the framework $\bar{ξ}$ is generated. Check whether $F (\bar{ξ}) < F (ξ_{min})$ , and update $ξ_{min} = \bar{ξ}$ . Otherwise, generate a stochastic variable $η$ that is uniformly distributed in the interval $(0, 1)$ . Calculate the value of $P (\bar{ξ} | ξ_{min}, T)$ which is the acceptance probability of $\bar{ξ}$ under the current temperature T and the optimal tentative point $ξ_{min}$

P (\bar{ξ} | ξ_{min}, T) = min {1, \exp ((F (\bar{ξ}) - F (ξ_{min})) / ζ T)}

(31)

If $η \leq P (\bar{ξ} | ξ_{min}, T)$ , update $ξ_{min} = \bar{ξ}$ . Otherwise, go to step 5.

Step 5: if $j < j_{max}$ , update $j = j + 1$ and go to step 2, otherwise, go to step 6.

Step 6: if the iteration termination condition, which could be controlled through the annealing temperature, is satisfied, the algorithm stops. Otherwise, update the temperature $T = ϑ T$ , where $0 < ϑ < 1$ is the cooling coefficient. Update $j = 0$ and go to step 2.

Numerical example

To evaluate the effectiveness of the proposed method, numerical experiments are conducted on a corridor (shown in Figure 6). This corridor could also be viewed as a simple railway network since the train routes on this corridor are interrelated. There are totally six different train routes and there are both high-speed trains and medium-speed trains on each route. The hollow vertices represent the stations at which a high-speed train is scheduled to stop. The solid vertices represent the stations at which a medium-speed train is scheduled to stop. The shaded vertices represent the stations at which both high-speed train and medium-speed train are scheduled to stop.

Figure 6.

Train routes in the basic test case.

Table 1 shows the data of the train operation scheme used in the basic experiment. $| K' |$ and $| K'' |$ are the number of high-speed trains and the number of medium-speed trains to be scheduled in each route, respectively.

Table 1.

Train operation scheme for the basic test case.

Route	First station	Last station	Mileage (km)	$\| K' \|$	Speed per hour (km/h, $k \in K'$ )	$\| K'' \|$	Speed per hour (km/h, $k \in K''$ )
1	SHHQ	BJN	1435	32	300	10	250
2	SHHQ	BBN	535	19	300	6	250
3	QFD	BJN	580	19	300	6	250
4	BJN	SHHQ	1435	32	300	10	250
5	BBN	SHHQ	535	19	300	6	250
6	BJN	QFD	580	19	300	6	250

Stations BJN, QFD, BBN, and SHHQ are able to fulfill the task of first-grade EMU maintenance. Other basic parameters for the construction of the EMU circulation plan are as follows: $α = 2880 min$ , $β = 4000 km$ , $ε = 15 min$ ,²⁵ and $ε = 120 min$ .²⁶ The minimum headway between consecutive arrivals and departures of trains at each station is set to be 4 min. The dwell time at each station is set to be 2 min. The available time range for trains to depart from the first station of their route is from 7:00 to 24:00.

When the total extra dwelling time of medium-speed trains generated due to being overtaken by high-speed trains and the total extra overtaking are viewed as penalties of the objective function of the main model, and the results of the initial solution and optimal solution obtained in the solving process of the proposed method are listed in Table 2. I-NLDT (O-NLDT) means the number of trains with long dwelling time in each route in the initial solution (optimal solution). I-R-NLDT (O-R-NLDT) means the ratio of I-NLDT (O-NLDT). I-NEO (O-NEO) means the number of trains with extra overtaking in the initial solution (optimal solution). I-R-NEO (O-R-NEO) means the ratio of I-NEO (O-NEO). I-HD means the initial solution’s homogeneity of the departure times of trains at the first stations. O-HD means the optimal solution’s homogeneity of the departure times of trains at the first stations.

Table 2.

Results of the initial solution and optimal solution in the basic case.

Route	I-NLDT	O-NLDT	I-R-NLDT	O-R-NLDT	I-NEO	O-NEO	I-R-NEO	O-R-NEO	I-HD (×10⁴)	O-HD (×10⁴)
1	3	3	0.3	0.3	2	0	0.2	0.0	–	–
2	1	1	0.17	0.17	0	1	0	0.17	40.2	27.9
3	3	1	0.5	0.17	1	0	0.17	0.0	6513.0	6322.9
4	2	1	0.2	0.1	3	1	0.3	0.1	–	–
5	0	0	0.0	0.0	1	0	0.17	0.0	7745.4	7485.6
6	2	1	0.33	0.17	3	1	0.5	0.17	43.5	29.3

It can be seen from Table 2 that the proposed method is able to filter the solutions with excessive overtaking or dwelling time due to being overtaken. The iteration process of the proposed method is shown in Figure 7. The total connection time reduces from 6.3859 × 10⁴ min of the initial solution to 5.7345 × ×10⁴ min of the optimal solution. The total travel time reduces from 4.2701 × 10⁴ min of the initial solution to 4.2015 × 10⁴ min of the optimal solution. The corresponding number of EMUs needed reduces from 74 of the initial solution to 69 of the optimal solution.

Figure 7.

Iteration processes of different solution methods for the basic case.

Figure 7 also makes comparison between the solving process of the proposed simulated annealing with managed neighborhood of first departure times (MD-AS) and the solving process of the general simulated annealing (AS) in which the neighborhood is unlimited. Obviously, the neighborhood of MD-AS is more reasonable than that of the AS.

Figure 8 displays the results of the branch-and-bound method obtained from some intermediate solutions in the process of solving the main model when the number of branches is set to be 4, 5, and 6. It can be seen from the figure that there are obvious differences between the results when the number of branches varies. Therefore, the searching space is broadened in the improved branch-and-bound method.

Figure 8.

Results of the branch-and-bound method when the branching number varies.

It is worth noting that the scale of the test case (e.g. the number of trains or the mileage of the railway) in our article is close to the real-world situations. The running time of the proposed approach in the experiment is about several hours. Specifically, the path selection at each iteration of the main model is within a few seconds. The construction of the EMU circulation plan for each iterative timetable is within a few minutes. Since the timetable and EMU circulation plan are obtained at a strategic level, the computation time of the proposed approach is reasonable for implementation in the practical applications.

Figure 9 displays the total connecting time of the intermediate points in the process of solving the second experiment case. The second case is constructed by reducing the number of trains in the basic case. Other parameters remain the same. Three types of results are shown for each intermediate point. The first type is like the tight connection. It can be seen from this type of results that when the connection relationships between trains are constructed by some heuristic rules, then the total connection time could be too long. The second type is the previous branch-and-bound method. The third type is the improved branch-and-bound method. It can be seen from the figure that the improved branch-and-bound method has obvious potential to further enhance the quality of the optimal circulation plan for each intermediate point. This feature of the improved branch-and-bound method plays an important role in helping the proposed framework find a balanced optimal solution.

Figure 9.

Comparison between the connecting times for intermediate solutions.

Other additional cases are constructed by reducing different number of high-speed trains and medium-speed trains in the basic case. Figure 10 makes the comparison between the total connection times of the random initial solution, the optimal solution obtained by the proposed method, and two optimal solutions obtained by the ordinary method. The ordinary method indicates that the optimal timetable and EMU circulation plan are found in a separated way. Specifically, the timetable with optimal efficiency is first obtained, and then the optimal EMU circulation plan is calculated based on this optimal timetable. Figure 11 makes the comparison between the total travel times of the random initial solution, the optimal solution obtained by the proposed method, and two optimal solutions obtained by the ordinary method. The connection times of the two solutions obtained by the ordinary method are calculated based on two optimal timetables whose efficiencies are nearly the same. It can be seen from Figures 10 and 11 that for each additional case, even if the total travel times of different timetables are close to each other, the connection times calculated based on them could have great difference. Figure 12 makes the comparison between the total objective values of the solutions of each case.

Figure 10.

Comparison between total connection times for different cases.

Figure 11.

Comparison between total travel times for different cases.

Figure 12.

Comparison between total objective values for different cases.

As seen from Figure 12, in the tests conducted for different cases, the total objective values of the output solutions of the proposed approach are on average 5.32% better than the solutions of the ordinary method and 8.87% better than the randomly generated initial solution. It indicates that the proposed method has a positive influence on the coordinate optimization of the timetable and the EMU circulation planning than the method that handles the problems in a sequential way.

Table 3 lists the total circulation times (TCT) of the randomly generated initial solution, the optimal solution obtained by the proposed framework for coordinate optimization of the timetable and EMU circulation plan (FCO), and two solutions obtained by the ordinary method (OM) for each case. According to TCT and equation (18), the number of EMU (NEMU) needed to fulfill all the train tasks for each case could be calculated. Table 3 also makes comparison between the NEMUs of different solutions. According to the results, the NEMUs of the solutions of the proposed approach is up to five (at least three) less than the initial solutions and up to four (at least one) less than the solutions of the method that handles the problems in a sequential way. It indicates that the proposed method is able to help the railway system save the EMUs which are expensive equipment.

Table 3.

Comparison between TCTs and NEMUs.

Case	Initial solution (×10²)		Solution-FCO (×10²)		Solution 1-OM (×10²)		Solution 2-OM (×10²)
Case	TCT	NEMU	TCT	NEMU	TCT	NEMU	TCT	NEMU
1	1008	70	964.8	67	993.6	69	993.6	69
2	1008	70	950.4	66	979.2	68	993.6	69
3	993.6	69	950.4	66	964.8	67	964.8	67
4	979.2	68	921.6	64	964.8	67	950.4	66
5	979.2	68	921.6	64	950.4	66	936	65
6	964.8	67	907.2	63	950.4	66	921.6	64
7	936	65	878.4	61	907.2	63	907.2	63
8	936	65	878.4	61	907.2	63	907.2	63
9	936	65	864	60	907.2	63	921.6	64

TCT: total circulation times; NEMU: number of electric motor train unit.

Table 4 lists the connection times of the DCM1-T, the initial feasible circulation plan (IFC) the optimal circulation plans obtained by branch-and-bound method, and the improved branch-and-bound method corresponding to different solutions of each case. INS represents the initial solution. It can be further seen from the results that the improved branch-and-bound method has stronger strength in finding an optimal circulation plan that is equal to or much close to the connection time of DCM1-T than the basic branch-and-bound method. Specifically, the connection times of solutions obtained by the proposed method is on average 1.13% more than the corresponding DCM1-Ts, on average 22.24% less than the solutions obtained by a tight-connection-based heuristic method, and on average 2.73% less than the solutions obtained by the original branch-and-bound method. It is worth mentioning that since DCM1-T does not consider the maintenance of the EMUs, the connection time of DCM1-T provides an approximate lower bound of the circulation planning problem as the submodel has described.

Table 4.

Comparison between the connection times corresponding to different solutions of each case.

Solution	Case	DCM1-T	IFC	B&B	M-B&B	Case	DCM1-T	IFC	B&B	M-B&B
INS	1	59196	83676	59196	59196	6	55659	72939	59979	57099
FCO		57137	68657	55697	55697		51751	67591	51751	51751
OM-1		58610	70130	58610	58610		53212	71932	58972	56092
OM-2		58540	75820	64300	58540		51911	70631	54791	53351
INS	2	57791	75071	63551	59231	7	55244	72524	55244	55244
FCO		52855	70135	60055	54295		49770	67050	49770	49770
OM-1		57294	80334	57294	57294		52702	67102	52702	52702
OM-2		57297	77457	64497	58737		52895	71615	52895	52895
INS	3	58272	75552	58272	58272	8	55350	72630	55350	55350
FCO		54751	72031	54751	54751		50037	70197	55797	50037
OM-1		56265	70665	56265	56265		52979	61619	52979	52979
OM-2		56240	69200	56240	56240		53013	68853	53013	53013
INS	4	57627	70587	57627	57627	9	51203	71363	55523	55523
FCO		51290	67130	55610	52730		48878	61838	48878	48878
OM-1		57097	72937	57097	57097		50540	64940	56300	53420
OM-2		55635	70035	61395	55635		54871	67831	54871	54871
INS	5	56604	76764	60924	58044
FCO		52969	70249	55849	52969
OM-1		54528	68928	61728	55968
OM-2		53014	70294	55894	54454

DCM1-T: Degree Constraint Minimum 1-Tree; IFC: initial feasible circulation plan; B&B: branch-and-bound method; M-B&B: improved branch-and-bound method.

Conclusion

In this article, a hierarchical framework is proposed for the coordinate optimization of train timetable and EMU circulation plan. The framework integrates the problem of timetabling and EMU circulation planning and consists of a main model and a submodel. The timetabling main model aims at finding an optimal timetable with balanced efficiency and the quality of the corresponding EMU circulation plan. For each iteratively generated timetable, the efficiency is computed directly in the main model and the quality of the EMU circulation plan is evaluated through the submodel.

A simulated annealing–based method which is embedded with an improved branch-and-bound algorithm is proposed for solving the problem. As the quality of the optimal solution of the proposed framework is not only related to the performance of the algorithm of the main model but also related to the performance of the algorithm of the submodel, each algorithm is designed effectively. The neighborhood system of the simulated annealing process is predetermined when solving the main model. This strategy helps the algorithm ignore the obviously unreasonable region when finding an intermediate solution. The improved branch-and-bound method enhances the capability of searching the optimal circulation plan in a wider space.

The numerical experiments show that for different test cases, the proposed method doses well in finding a timetable which not only has the optimal efficiency but also provides a good basis for the constructing of an optimal circulation plan. Furthermore, the improved branch-and-bound method is able to ensure that the circulation plan of each intermediate timetable reaches its approximate lower bound.

In our future research, we will extend the proposed framework that incorporates the robustness into the output solutions. How to preserve the stability of the timetable and the EMU circulation plan against the stochastic disturbances (e.g. Markov-type disturbances) needs to be further studied.

Footnotes

Academic Editor: Hamid Reza Karimi

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was financially supported by the Science and Technology Research and Development Program of China Railway Corporation (KTD16011531).

References

Caprara

Fischetti

Toth

Modeling and solving the train timetabling problem. Oper Res 2002; 50: 851–861.

Brännlund

Lindberg

Nõu

. Railway timetabling using Lagrangian relaxation. Transport Sci 1998; 32: 358–369.

Zhou

Zhong

Bicriteria train scheduling for high-speed passenger railroad planning applications. Eur J Oper Res 2005; 167: 752–771.

Lee

Chen

A heuristic for the train pathing and timetabling problem. Transport Res B: Meth 2009; 43: 837–851.

Castillo

Gallego

Ureña

. Timetabling optimization of a mixed double- and single-tracked railway network. Appl Math Model 2011; 35: 859–878.

Cacchiani

Furini

Kidd

MP.

Approaches to a real-world train timetabling problem in a railway node. Omega 2016; 2016: 97–110.

Chu

Zhang

Chen

. Pareto optimal train scheduling for urban rail transit using generalized particle swarm optimization. Adv Mech Eng 2016; 8: 1–15.

Zhao

Yang

AZ.

Research on usage of high speed passenger trains on uncertain railroad region. J China Railw Soc 1997; 19: 15–19.

Peeters

Kroon

Circulation of railway rolling stock: a branch-and-price approach. Comput Oper Res 2008; 35: 538–556.

10.

Alfieri

Groot

Kroon

. Efficient circulation of railway rolling stock. Transport Sci 2006; 40: 378–391.

11.

Cacchiani

Caprara

Toth

A Lagrangian heuristic for a train-unit assignment problem. Discrete Appl Math 2013; 161: 1707–1718.

12.

Abbink

Bianca

VDB

Kroon

. Allocation of railway rolling stock for passenger trains. Transport Sci 2004; 38: 33–41.

13.

Shi

Zhou

Wei

. Optimized model and algorithm of motor train-sets scheduling for dedicated passenger lines. J China Railw Soc 2011; 33: 8–13.

14.

Zhao

Tomii

An algorithm for train-set scheduling on weekday based on probabilistic local search. Syst Eng 2004; 24: 123–128.

15.

Cacchiani

Toth

Nominal and robust train timetabling problems. Eur J Oper Res 2012; 219: 727–737.

16.

Cacchiani

Caprara

Galli

. Railway rolling stock planning: robustness against large disruptions. Transport Sci 2012; 46: 217–232.

17.

Wei

Qiu

Karimi

HR.

Reliable output feedback control of discrete-time fuzzy affine systems with actuator faults. IEEE T Circuits: I 2017; 64: 170–181.

18.

Wei

Qiu

Karimi

. Model reduction for continuous-time Markovian jump systems with incomplete statistics of mode information. Int J Syst Sci 2014; 45: 1496–1507.

19.

Wei

Qiu

Karimi

. Filtering design for two-dimensional Markovian jump systems with state-delays and deficient mode information. Inform Sciences 2014; 269: 316–331.

20.

Wei

Qiu

Karimi

. New results on H_∞ dynamic output feedback control for Markovian jump systems with time-varying delay and defective mode information. Optim Contr Appl Met 2014; 35: 656–675.

21.

Wei

Qiu

Karimi

HR.

Quantized H_∞ filtering for continuous-time Markovian jump systems with deficient mode information. Asian J Control 2014; 17: 1914–1923.

22.

Fisher

ML.

A polynomial algorithm for the degree-constrained minimum K-tree problem. Oper Res 1994; 42: 775–779.

23.

Zhou

Yue

. A branch and bound algorithm for the exact solution of the problem of EMU circulation scheduling in railway network. Math Probl Eng 2016; 2016: 1–14.

24.

Alexandre

SDC

Abilio

. Algorithms for the degree-constrained minimum spanning tree problem. Electron Notes Discret Math 2005; 19: 403–409.

25.

Han

Zhang

. Research on optimization model and algorithm of EMU circulation plan. J China Railw Soc 2013; 35: 1–8.

26.

Miao

Wang

Yang

ZX.

Research on the optimization of EMU circulation based on optimized connecting network. J China Railw Soc 2010; 32: 1–7.